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finite  element  analysis.  The  concepts  used  in  developing  the 
variational  statement  are  somewhat  different  from  those  used  by 
most  other  Investigators  and  appear  to  offer  certain  advantages 
for  inelastic  formulations.  Finally  results  obtained  with  the 
finite  element  analysis  are  compared  to  known  solutions  with 
good  results. 

For  the  convenience  of  the  reader  the  total  report  on  the  j 
project  is  presented  in  four  parts.  As  noted  above  a  description 
of  the  consolidation  theory  and  certain  theoretical  features  of 
the  finite  element  analysis  are  described  in  the  body  of  the 
main  report  (CR  84.006).  The  second  part  (CR  84.007)  describes 
the  numerical  evaluation  of  the  incremental  form  of  the  bounding 
surface  model.  Finally ^user's  manuals”  for  the  2-D  and  3-D 
finite  element  programs  are  given  in  two  additional  reports 
(CR  84.008  and  CR  84.009). 
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INTRODUCTION 


During  the  past  four  years,  extensive  theoretical,  experimental  and 
numerical  work  has  been  done  at  UCD  on  a  bounding  surface  plasticity  model 
for  cohesive  soils.  The  work  has  been  conducted  in  part  under  the  sponsorship 
of  contracts  and  grants  from  the  U.S.  Army  Waterways  Experiment  Station  (Re: 
Contract  #DACA  39-79-M-0059),  the  Civil  Engineering  Laboratory,  Naval 
Construction  Battalion  Center  (Re:  Orders  N62583-80-M-R478,  N62583-81-M- 
R320,  N62474-S2-C-8276  and  N62583-83-M-T062),  and  NSF  (Re:  Grants 
NSF-CME-79- 10835  and  NSF -CEE-82- 16995).  To  date  this  work  has  resulted  in 
the  reports  and  papers  listed  in  references  [1  -*■  14]. 

Although  continued  research  will  most  certainly  bring  refinements  and 
new  areas  of  application,  the  model  has  now  reached  a  stage  of  development 
where  it  can  be  used  as  a  practical  engineering  tool.  Such  use  requires  that 
it  be  incorporated  into  new  and  existing  finite  element  codes  for  the  analysis 
of  earth  structures.  In  order  to  simply  and  inexpensively  achieve  such  numerical 
implementations,  the  model  has  been  coded  into  a  master  subroutine  CLAY  and 
several  supporting  sii>routines  [2,5]. 

Because  of  the  evolutionary  nature  of  the  development  of  this  code,  with 
time  it  has  become  increasingly  inefficient  and  unwidely.  Thus,  to  improve  its 
readability,  to  incorporate  some  recent  minor  revisions  in  the  model  and  its 
numerical  implementation,  and  to  take  advantage  of  the  structured  nature  of 
FORTRAN  77,  the  subroutines  have  been  recoded.  In  this  process  it  was  found 
to  be  convenient  to  slightly  modify  the  calling  instructions  for  the  subroutine 
and  thus  there  is  a  need  to  revise  the  instructions  for  its  implementation  as 
given  in  [5]. 

It  is  the  purpose  of  this  report  to  discuss  the  recent  minor  changes  in 
the  model,  and  its  numerical  implementation,  and  to  give  detailed  instructions 


for  its  incorporation  into  new  and  existing  finite  element  codes.  In  addition, 
the  program  (EVAL)  [10]  written  to  evaluate  homogeneous  test  results  for 
cohesive  soils  is  discussed.  This  report  is  for  the  most  part  a  replacement  for 
Section  2  and  Appendices  I  and  n  of  f 51.  Frequent  reference  will  be  made  to 
equations  and  figures  from  [5];  such  references  are  expressed  in  the  form  (5-N), 
etc. 

1.  Recent  Modifications  to  the  Bounding  Surface  Plasticity  Model  for  Cohesive 
Soils 

This  report  is  not  intended  to  stand  alone,  but  rather  to  be  a  replacement 
for  certain  sections  of  [5],  Thus,  for  a  comprehensive  description  of  the 
underlying  theory,  notation  and  the  background  for  this  report,  the  reader  is 
referred  to  [5]. 

In  this  section  the  modifications  that  have  been  made  to  the  bounding 
surface  plasticity  model  for  cohesive  soils  since  the  publication  of  [5]  are 
discussed. 

Modification  1: 

Following  a  recent  paper  by  Ha  fa  lias  [15],  the  "radial"  mapping  rule  of 
eq.  (5-12)  has  been  generalized  to  give  an  image  stress  state  of 

T  =  B(I  -  Ic)  ♦  Ic,  Sj.  =  8  s..,  3  =  BJ,  a  =  a  (1) 

This  generalization  has  introduced  the  possibility  of  using  a  projection  center 
(lc  =  CI0)  in  stress  space  different  than  the  origin.  This  is  illustrated  in  Fig. 
1;  note  that  the  parameter  "C"  has  no  relation  to  point  "C"  of  Figure  5-2. 
For  C  =  0  one  retrieves  the  previous  formulation.  This  modification  introduces 
a  kind  of  "hydrostatic  back  stress,"  Ic,  and  allows  for  the  prediction  of  immediate 
negative  pore  pressure  development  for  heavily  overconsolidated  samples  under 
undrained  loading  conditions.  With  the  projection  center  at  the  stress  origin, 
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one  would  always  have  an  initial  positive  porewater  pressure  built-up  even  for 
large  OCR  values. 

Modification  2: 

The  denominator  of  F -  6  in  eq.  (5-53)  has  been  replaced  by  the  quantity 
<r  -  s6>  with  s  being  an  elastic  nucleus  parameter  and  the  symbol  F  being  now 
simply  written  as  r.  Thus,  whenever  6  >  r/s  the  brackets  yield  a  zero  value 
and  according  to  eq.  (5-28)  Kp  =  «.  Observe  that  Kp  <»  as  <5  -*■  r/s.  A 
geometric  interpretation  of  this  behavior  is  that  s  indirectly  defines  a  stress 
domain  of  purely  elastic  response  within  the  bounding  surface;  it  is  not  necessary 
to  interpret  the  boundary  of  this  zone  as  a  yield  surface  (no  associated  consistancy 
condition,  etc.).  This  domain  is  called  the  ’’elastic  nucleus"  and  is  shown  in 
Fig.  1.  The  quantity  s  which  controls  the  size  of  the  elastic  nucleus  can  be  a 
function  of  the  state  including  the  accumulated  deviatoric  plastic  strains.  The 
elastic  nucleus  controls  the  possible  stabilization  of  cyclic  undrained  stress  paths, 
allowing  for  stabilization  or  failure  depending  upon  the  amplitude  of  the  cyclic 
deviatoric  stress. 

Modification  3: 

Eq.  (5-55)  is  replaced  by 


<1. 


V  +  h 


X  —  tc 


(2) 


This  expression  predicts  that  for  tensile  stresses  the  soil  looses  all  cohesion  at 
a  maximum  but  finite  void  ratio.  The  original  expression,  eq.  (5-55),  implied  a 
plastic  void  ratio  tending  towards  infinity  as  I  0.  (Recall  that  IQ  is  the 
internal  variable  controlling  the  size  of  the  bounding  surface  and,  in  general,  is 
not  the  current  value  of  1). 

Modification  »; 

The  quantity  pa  (atmospheric  pressure)  in  eq.  (5-84)  has  been  replaced  by 
(1  *  “  ^l>  +  term  l  is  now  used  (instead  of  p  )  to 
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non-dimensionalize  h(a)  and  the  term  (1  +  e  )/(A  -  <)  is  introduced  for  reasons 

o 

of  convenience  since  Kp  depends  on  this  factor  as  well.  The  use  of  lQ  in  this 
expression  renders  all  undrained  stress  paths  at  a  given  OCR  similar  when 
normalized  with  respect  to  1Q.  If  natural  strains  are  used  instead  of  engineering 
strains,  the  inital  void  ratio  (eQ)  is  replaced  by  the  current  value  (e). 
Modification  5; 

The  term  h(a)  in  eq.  (5-84)  has  been  replaced  by 

h(a,z)  =  zf"  hj(a)  +  (l  -  zJ")  h? 
where  (see  Fig.  1  and  refer  to  eq.  (5-58)) 

h,(ot)  -  - : — 5 —  h 

1  1  +  y  -  (1-y)  sm3a  c 

h 

W  s  ,  with  he,  hc  being  the  values  of  h  for  triaxial  compression  and  extension 
c 

respectively. 

..3  3  JR 
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m  =  a  constant  which  applies  to  both  extension  and  compression  and  is  thus 
not  a  function  of  the  Lode  angle  a. 

hj  =  the  shape  hardening  parameter  for  states  on  the  I-axis  (i.e.  for  z  =  0). 

This  term  assures  continuity  when  crossing  the  I-axis;  the  modification  was 
incorporated  to  improve  numerical  behavior  in  this  region,  however,  it  has  very 
little  affect  on  the  model  predictions. 

Incorporating  modifications  2,  4  and  5,  eq.  (5-84)  becomes 

KP  -  %  ♦  xrr  [<Io  -  V  ♦  hi  +  (1  - 

i  &  (3) 

f9  f,t  ♦  } 
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The  model  as  described  in  [5]  and  subject  to  the  above  modifications  gives 
a  comprehensive  description  of  the  three-dimensional  stress-strain  behavior  of 
cohesive  soils.  The  practical  use  of  the  model  depends  upon  its  numerical 
implementation  in  new  and  existing  finite  element  codes  for  earth  structures; 
this  topic  will  be  discussed  in  the  following  section. 


2.  NUMERICAL  IMPLEMENTATION* 

2.1  fncrementalization  of  the  Bounding  Surface  Plasticity  Rate  Equations 

The  bounding  surface  plasticity  theory  for  cohesive  soils  is  expressed  in 
terms  of  effective  stress,  whereas  most  soil  related  problems  involve  the 
application  and  calculation  of  total  stress.  The  difference  between  the  total 
and  effective  stress  is  simply  the  pore  water  pressure  u.  This  section  deals 
with  the  incremental  relation  between  the  strain  and  effective  stress  components; 
the  pore  water  pressure,  and  hence  the  total  stress,  is  dealt  with  in  a  following 
section. 

Factoring  the  strain  increment  from  eq.  (5-30)  gives: 


°ij  =  Dijk*  ek l 
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This  section  serves  as  a  replacement  for  the  corresponding  section  in  [51. 


where 


L  >  0 
L<  0 


( load ing) 
(unloading) 


(6) 


L  = 


l 

D 


G.  . 

I  J 


+  /Tg  f,  f/s*  Ski 

3  cos  3a  ’a  \  n2 


n  ■  * 9  K  (f’t J  * G  (F'J*  I  N' 


(7) 

(8) 


Eq.  (4)  relates  the  tensor  components  of  (effective)  stress  and  strain.  For 
finite  element  analysis  purposes,  it  is  more  convenient  to  express  this  relationship 
in  matrix  form,  i.e., 

{a}  =  fr>]  {1}  (9) 

T 

where  ([  1  is  the  matrix  transpose) 


MT  *  (ax,  ayt  az,  txy,  tx2,  xyz)  (10) 

{el  =  <ex»  V  V  \y>  ^Xz’  V 

The  tensor  components  of  shear  strain  are  one-half  of  the  engineering 
components  y...  The  [Dl  matrix  is  related  to  the  components  of  the  D.M  „ 

IJ  IJKX- 

tensor  as  follows  (because  of  the  symmetry  of  the  stress  and  strain  tensors, 


interchanging  i  and  j  or  k  and  2.  in  eq.  (5)  results  in  the  same  quantity): 
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In  order  to  use  eq.  (9)  in  a  finite  element  program,  it  must  be  expressed 
in  an  incremental  form.  Consider  the  Nth  step  of  an  incremental  analysis;  i.e., 
the  solution  has  been  found  at  N-l,  and  it  is  now  desired  to  calculate  the 
incremental  change  that  will  give  the  solution  at  N. 

Numerous  numerical  methods  have  been  developed  for  calculating 
incremental  solutions  to  plasticity  problems.  A  simple  family  of  solutions  which 
can  be  viewed  as  approximations  to  Newton-Raphson's  method  will  be  discussed 
here  [12,16];  members  of  this  family  include  the  classical  method  of  "successive 
approximations"  and  the  popular  "tangent  stiffness  method"  [16].  With  little 
additional  effort  the  equations  to  be  discussed  can  be  adapted  for  use  with  a 
number  of  other  methods. 

Because  of  the  nonlinear  nature  of  elastic-plastic  behavior,  iteration  is 
in  general  required  to  establish  the  incremental  change.  At  the  end  of  the  K-l 
iteration,  the  estimates  of  the  stress  and  strain  states  at  N  are  given  by  the 
expressions: 

Mn,k-i  *  (“Wi  *  !toW-t  un 

Ic>n,k-i  =  ,e,n-i  *  H,k-i  <12) 

The  iteration  process  is  continued  until  some  specified  convergence 
criterion  is  satisfied,  or  until  a  specified  maximum  number  of  iterations  have 
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failed  to  yield  convergence  (an  indication  of  possible  failure  and/or  unstable 
behavior). 


Even  though  rate  independent  behavior  is  being  considered,  it  is  convenient 
to  think  in  terms  of  the  time  history  of  the  quantities  involved.  If  eq.  (9)  is 
integrated  from  time  t^j  to  t^,  it  yields: 

*N  tN 

/  {0}  dt  =  r  ml  {e}  dt  (13) 

Vi  lN- 1 

or  {Ao},^  =  jT  P3l  {c}  dt  (14) 

Vl 

A  simple  two  point  formula  is  used  to  approximate  the  above  integral, 

(0  <_  0,  <_  1): 

|to)N  =[u  -  V  ©1N.,  *  e,  [d!n]  lsclN  05) 

Values  of  0j  =0,  1/2  and  1  give,  respectively,  forward,  trapezoidal  and  backward 
integration.  Although  trapezoidal  integration  is  most  accurate,  in  rare  instance 
it  may  be  advantageous  to  use  9  t  =  0.0  in  order  to  reduce  iteration  requirements 
(at  the  expense  of  smaller  step  sizes). 

Because  [OJN  and  [nJN  j  are  functions  of  the  stress  and  strain  states 
at  N  (see  eqs.  (5,6)),  it  is  necessary  to  base  their  values  on  the  stress  and  strain 
estimates  of  the  previous  iteration  (see  eqs.  (11)  and  (12)).  The  fact  that 
is  possibly  a  function  of  { Ao}N  and  {Ae)N  requires  some  explanation.  The 
dependence  is  a  result  of  the  discontinuous  nature  of  plasticity  behavior  at  the 
initiation  of  an  unloading  process  (see  eq.  (6)).  For  the  simple  one-dimensional 
example  shown  in  Figure  2,  the  stiffness  will  be  O'  or  D"  depending  upon 

whether  the  step  is  to  N*  or  N".  While  in  a  one-dimensional  problem  this  would 
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usually  be  known  a-priori,  for  a  finite  element  that  is  a  part  of  a  complicated, 
highly  statically  indeterminate,  two  or  three-dimensional  structure  it  can  only 
be  established  by  the  iteration  process.  If  this  dependence  of  lDlN_j  on  the 
iteration  process  is  ignored  and  6  j  is  taken  as  zero,  then  iteration  is  not 
necessary,  however,  very  small  time  steps  are  required  for  accuracy;  in  general, 
this  procedure  is  not  recommended.  The  predicted  value  for  [D]^  is  denoted 
by  lD]|yj£  j  etc.  The  equation  resulting  from  substituting  these  estimates  into 
eq.  (15)  is  used  to  relate  the  estimates  of  { Ao}^  and  {Ae}^  for  iteration  K, 
i.eM 

^N,K  =  ^N,K-1  ^N,K  (16) 

where 

^N,K-1  =  t(1  '  0l)CDlN-l,K-l  +  61  (l7) 

Cq.  (16)  is  the  desired  incremental  stress-strain  relation  for  iteration  K 
of  increment  N,  and  in  the  Newton-Raphson  method,  is  used  for  the  calculation 
of  the  residual  vector  [121.  In  addition  the  Newton-Raphson  method  requires 
the  Jacobian  matrix.  An  approximation  to  the  Jacobian  can  be  written  in  the 
form  [12,14]. 

"  e2)lDlN-l,K-l  +  02t^IN,K-l1  (I8) 

Values  of  6  j  of  1/2  and  1  correspond  respectively  to  the  methods  of  successive 
approximation  and  tangent  stiffness  [12,16].  The  simple  test  evaluation  program 
EVAL  discussed  in  [10]  and  a  later  section  of  this  report  uses  62  =  1/2,  the 
consolidation  codes  discussed  in  [14]  permit  02  to  be  specified  by  the  user;  all 
these  applications  use  8|  =  1/2  in  eq.  (17). 

Thus,  for  finite  element  implementation  what  is  required  is  the  calculation 
of  the  matrices  and  Th«  combining  of  these  matrices 

in  eq.  (17)  and  (IS)  (or  their  use  in  some  other  fashion  for  a  different  nonlinear 


solution  algorithm)  is  left  for  the  finite  element  program  which  calls  the  master 
$ti>routine  CLAY  that  has  been  written  for  their  evaluation.  The  main  goal  of 
this  report  is  a  description  of  the  steps  necessary  for  incorporating  this  subroutine 
into  new  or  existing  finite  element  codes  (Section  2.3).  The  calculation  of 
[n]N  j  ^  j  and  [F>  1  |  proceeds  directly  from  eq.  (5)  with  only  a  few  steps 

needing  elaboration. 

The  parameter  9  in  eq.  (5-50)  is  undefined  for  a  zero  value  of  the  first 

effective  stress  invariant  I.  The  numerical  problems  associated  with  a  zero  or 

near  zero  value  of  I  are  avoided  by  arbitrarily  replacing  the  value  of  |I|  by 

lCT^P  for  such  cases  (where  P  is  atmospheric  pressure);  the  corresponding 
a  a 

expression  in  the  equation  of  "Modification  5"  (Section  1)  is  treated  in  a  similar 
fashion.  In  general,  this  arbitrary  action  does  not  significantly  affect  the 
calculated  properties. 

As  the  soil  state  approaches  the  bounding  surface,  a  stress  state  outside 
of  the  surface  may  be  predicted  in  a  particular  iteration.  Because  such  a 
prediction  has  no  meaning,  the  state  is  assumed  instead  to  fall  on  the  surface; 

r 

i.e.,  8  is  restricted  to  be  >  1.0  and  — - y-  is  restricted  to  be  >  0. 

At  any  instant,  the  size  of  the  bounding  surface  is  determined  by  the 
value  of  1  (see  Figure  1).  The  differential  change  in  I  is  obtained  from  eq. 
(2),  i.e., 

d,0  -  nir  (<lo  -  V  *  '«.'(dekk  •  k d,)  «’> 

Two  cases  must  now  be  considered  in  integrating  (19):  If  I  >  1^,  eq.  (19) 
becomes: 

d,0  ■  (1*'o)  "n!-  <d£kk  -  kdl) 


(20) 
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Dividing  by  I  and  integrating  the  resulting  expression  for  increment  N  gives: 
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The  first  two  integrals  may  be  evaluated  exactly,  while  the  third  is  approximated 
by  the  trapezoidal  rule  giving: 


in 
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If  I  <_  I^,  eq.  (19)  becomes 
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TTr"r'kkN,K-i'  *  <jw  *  Mw-' 


(22) 


(23) 


dl  s  (1+e  ) 


TZk  (d€kk  -  5R 


dl) 


(24) 


integrating  the  resulting  expression  for  increment  N  gives: 


Vk-1  ( 1+e  )  \  *?N,K- 

/  d'o  ■  T^-  *»  / 

Vi  K.i 


^,K-I 


de 


kk 


*N,K- 1 

/k 

Vi 


dl 


(25) 


Evaluating  the  first  two  integrals  exactly  and  approximating  the  third  by  the 
trapezoidal  rule  gives: 


(1+e  ) 

I  -  j  +  0  j 

°N,K-l  Vl  X~K  1 


4€kW  *  (vi  *  ^.k-!  4'n-k-1 


(26) 


Because  the  point  of  switching  from  eq.  (20)  to  eq.  (24)  (controlled  by  the  value 
of  I^)  is  somewhat  arbitrary,  consideration  was  not  given  to  the  possibility  of 
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this  changeover  occurring  in  mid  increment.  Instead  the  value  of  I  is  used 

°N-1 

to  make  the  decision  for  the  whole  increment. 

UntU  convergence  occurs,  the  estimates  of  { Ao}N  and  {  Ae}N  K  j  used 
in  the  calculation  of  the  incremental  properties  [D]n  K  1  do  not  in  fact  satisfy 
the  incremental  stress-strain  relation,  eq.  (16).  Because  this  inconsistency 
disappears  as  global  convergence  occurs  (i.e.,  as  [HlNK  j  a  [fi]N  K  2^  lt  is 
not  absolutely  necessary  to  take  any  special  steps  to  avoid  it,  however,  numerical 
experimentation  has  indicated  computational  advantages  in  doing  so  [12].  Thus, 
local  (i.e.,  within  subroutine  CLAY  each  time  it  is  called)  iteration  is  introduced 
in  the  calculation  of  [D]  to  remove  the  inconsistency  [17].  Using  {Ae}^K  j 
and  {  Ao|n  k_j*  n  calcu,ate<1*  “Hte  values  of  fAe)NK  j  and  [f)JNK  j 

are  then  used  in  conjunction  with  eq.  (16)  to  calculate  a  new  estimate  of  stress 
(Ao}^  K1  which  is  in  turn  used  with  {Ac)NK  l  to  calculate  a  new  estimate 
for  the  incremental  properties  This  process  is  continued  until 

convergence  is  achieved  for  Because  of  the  global  iteration  (the 

basic  iterative  procedure  of  the  calling  finite  element  program)  which  also  tends 
to  remove  the  inconsistency,  it  appears  that  the  convergence  limit  on  the  local 
iteration  can  be  considerably  less  restrictive  than  the  global  requirement.^  The 
stress  estimate  is  iteratively  modified  (instead  of  the  strain  estimate)  in  order 
to  maintain  a  compatible  global  displacement  field  as  required  by  the  admissibility 
conditions  of  the  finite  element  procedure.  The  introduction  of  local  iteration 
(for  ail  points  in  the  finite  element  grid  where  the  incremental  stress-strain 
properties  are  required)  of  course  substantially  increases,  in  a  given  global 
iteration,  the  computational  time  for  subroutine  CLAY,  however,  there  is  a 
corresponding  reduction  in  the  number  of  global  iterations. 


In  [14],  ten  times  the  global  limit  is  used. 
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2.2  Calculation  of  Pore  Voter  Pressure 

Many  finite  element  programs  have  special  provisions  for  handling  pore 
water  pressure  calculations  [14],  For  such  situations,  the  soil  properties 
sii>  routine  need  only  supply  a  relation  between  the  increments  of  strain  and 
effective  stress,  i*.,  eqs.  (16)  and  (18).  For  such  cases  the  remainder  of  this 
section  has  no  significance. 

There  are  three  possibilities  concerning  the  development  of  pore  water 
pressure  in  soil:  ideal  drained  conditions  (where  the  pore  water  pressure  is 
identically  zero),  ideal  undrained  conditions  (where  the  soil  is  completely 
saturated,  and  no  flow  of  water  occurs  whatsoever),  and  the  more  realistic 
situation  where  there  is  a  global  flow  of  water  and/or  the  filling  of  voids.  In 
many  analyses  ideal  drained  or  undrained  conditions  are  assumed,  even  though 
they  may  only  be  approximately  true. 

The  total  stress  rate  <?.  is  the  sum  of  the  effective  stress  rate  and  the 
pore  water  pressure  rate: 

o?.  =  o. .  +  u  6..  (27) 

U  «)  i) 

For  drained  conditions  u=0  and  o'.  =  o..,  and  eqs.  (16)  and  (18)  are  the 
desired  relations  between  the  total  stress  increment  and  the  strain  increment. 

For  undrained  conditions  there  are  several  possible  ways  of  proceeding. 
The  traditional  approach  has  been  to  neglect  the  (slight)  compressibility  of  the 
water  and  the  soil  particles,  and  thus  assume  incompressible  material  behavior. 
However,  the  finite  element  analysis  of  incompressible  materials  requires  a 
special  formulation  [18,19,20], 

In  order  to  avoid  having  to  deal  with  separate  formulations  for  drained 
and  undrained  conditions,  it  is  often  convenient  to  express  them  in  a  common 
form  (the  numerical  consequences  of  this  step  are  discussed  below).  This  can 


be  accomplished  if  the  slight  compressibility  of  the  soil  particles  and  the  pore 

* 

water  is  recognized  [21]  .  Thus,  the  pore  water  pressure  u  is  written  in  terms 
of  the  combined  bulk  modulus  r  of  the  soil  particles  and  the  pore  water  and 
the  resulting  (very  small)  volume  change  (note  that  as  T  *  •  the  soil 
becomes  incompressible,  and  that  drained  conditions  are  obtained  when  r=0): 


u  =  r 


For  undrained  conditions  the  value  of  T  is  very  large  compared  to  the 
terms  in  j.  Thus,  the  soil  behaves  as  a  "nearly  incompressible  solid" 

[18.19] ,  and,  consequently,  care  must  be  exercised  to  avoid  numerical  round-off 
and  element  "locking"  problems.  Two  approaches  are  commonly  used  to  achieve 
this  goal.  One  method  is  to  use  the  special  formulation  given  in  references 

[18.19]  for  incompressible  and  nearly  incompressible  solids,  while  the  other  is 
to  use  "reduced"  or  '•selective-reduced"  integration  [22,23]  for  the  element 
stiffness  matrix  (the  importance  of  selecting  a  proper  element  type  is  discussed 
in  [20,22,24]).  In  the  latter  case,  eq.  (28)  is  used  to  eliminate  u  from  eq.  (27); 


°ii  =0i)  *r 


Integration  over  increment  N  gives: 


=  *%,  * r  “S.  *>■ 

Using  the  above  equation  to  eliminate  Ac..  from  eq.  (16)  yields: 


lAot,N,K  s  ®*^N,K-1  tAelN,K 


An  alternative  interpretation  of  this  method  is  to  consider  the  undrained  soil 
as  incompressible,  and  to  incorporate  the  incompressibility  condition  by  means 
of  a  "penalty  function"  [22 1.  The  associated  "penalty  number"  corresponds 
to  the  bulk  modulus  I*. 
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where  (all  components  of  [d]  are  zero,  except  dn=<l22=c*33=^: 

'"'W-l  *  ®N,K-1  *  fdl  °2’ 

A  similar  procedure  is  followed  for  eq.  (18).  It  should  be  noted  that  eq.  (31) 
is  theoretically  valid  for  all  situations,  including  drained  conditions  (T=0). 
However,  for  very  large  values  of  T  (undrained  conditions),  problems  may  arise 
if  special  precautions  are  not  used  (201. 

2.3  Finite  Clement  Applications  -  Properties  Subroutine  (CLAY) 

A  properties  subroutine  (FORTRAN  77)  CLAY  (and  associated  subroutines) 
has  been  prepared  which  evaluates  the  incremental  stress-strain  properties  given 
by  the  bounding  surface  plasticity  model  for  cohesive  soils,  i.e.  the  matrices 
[0 1  jsj_i  k_i  and  which  appear  in  eqs.  (17)  and  (18).  A  listing  of  the 

subroutine  is  provided  in  Appendix  in.  The  subroutine  is  intended  for  incorporation 
into  new  or  existing  finite  element  programs  for  earth  structures;  such 
applications  by  the  authors  have  proven  to  be  straightforward  and  successful 
(12,141. 

For  each  iteration  of  each  increment,  and  for  all  points  (e.g.,  element 
centers  or  quadrature  points)  in  the  body  where  the  incremental  properties  are 
required,  the  parent  finite  element  program  calls  subroutine  CLAY.  The  call 
is  as  follows: 

CALL  CLAY  (IDIM,  INC,  ITNO,  ERMAX,  PROP,  STOR,  SIGB,  EPB,  DSIG,  HEP, 
DB,  DE,  UB,  DLTAU,  GAM,  KIND,  LARGE,  LOGIT,  THl). 

The  quantities  IDIM,  INC,  ITNO,  KINO,  LARGE  and  LOGIT  are  integer 
variables,  ERMAX,  UB,  DLTAIJ,  GAM,  and  THl  are  floating  point  variables,  and 
PROP,  STOR,  SIGB,  EPB,  DSIG,  DEP,  DB  and  DE  are  floating  point  arrays  of 
dimensions  (19),  (7),  (6),  (6),  (6),  (6),  (6,6)  and  (6,6),  respectively.  With  the 
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exception  of  IDIM,  which  is  discussed  in  a  later  paragraph,  the  arguments  in 
the  call  are  described  below: 


INC: 

ITNO: 

ERMAX: 

PROP: 


STOR: 


Increment  number  (the  first  increment  must  he  numbered  1) 
Iteration  number  (the  first  iteration  of  each  increment  must  be 
numbered  1) 

Convergence  limit  for  the  local  iteration,  suggested  values  are  .1  .01. 

An  array  containing  the  values  of  the  material  parameters  which 

describe  the  bounding  surface  plasticity  model  for  the  soil  at  the 

point  in  the  structure  for  which  the  incremental  properties  are  sought. 

The  parent  finite  element  program  must  read  and  store  the  values 

of  the  soil  parameters  for  each  different  type  of  soil  in  the  earth 

structure,  and  then,  for  each  call  to  CLAY,  present  the  appropriate 

values  for  the  element  in  question.  The  soil  parameters  are  stored 

in  the  array  in  the  following  order  (the  significance  of  the  various 

* 

parameters  are  described  in  detail  in  [5,131  ;  they  are  summarized 
in  Appendix  I):  X,  k,  Mc>  Rc,  Ac>  T,  P^,  v  or  G,  T,  m,  hc,  h2, 
n,  v,  r,  a,  c,  s).  That  is,  PROP(l)=X,  PROP(2)=ic,  etc.  At  the  time 
the  properties  are  read  into  the  main  program,  and  before  they  are 
stored,  and  must  be  multiplied  by  /3/9,  and  3,  respectively. 
It  is  suggested  that  subroutines  RPROP  and  TCHECK,  listed  with 
EVAL  in  Appendix  IV,  be  incorporated  into  the  parent  finite  element 
program  for  the  purpose  of  reading,  echo  printing  and  scaling  :he 
material  parameters.  The  input  formats  for  RPROP  are  those  given 
in  "III.  Material  Properties"  of  Appendix  I. 

This  array  is  used  to  store  certain  quantities  which  vary  with  the 

current  state  of  the  soil  (such  as  the  current  value  of  l  )  and  thus, 

o 


* 


See  beginning  of  the  report  for  a  discussion  of  modifications  to  the  model. 


IS 


are  for  a  given  step  in  the  analysis  unique  to  the  point  in  the  earth 

structure  under  consideration.  The  values  in  STOR  must  be  stored 

(after  each  call  to  CLAY)  by  the  parent  finite  element  program  for 

each  point  in  the  earth  structure  for  which  the  incremental  properties 

are  needed  (e.g.,  element  centers).  Prior  to  each  call  to  CLAY  the 

appropriate  values  for  the  point  in  question  are  retrieved  from  storage 

(i.e.,  from  a  two-dimensional  array  or  a  disk  file  which  stores  the 

values  for  each  element  in  the  system)  and  presented  to  the  subroutine. 

At  the  beginning  of  the  analysis  the  parent  finite  element  program 

must  initialize,  for  each  point  in  question,  STOR(l)  and  STOR(7),  with 

the  initial  value  of  3*P  and  e  respectively.  P  and  e  are  the 

o  o  o  o 

initial  values  of  the  preconsolidation  pressure  (Iq=3*Po  is  the  internal 
variable  controlling  the  size  of  the  bounding  surface,  Figure  1,  [2]) 
and  void  ratio.  These  values  are  not  read  by  subroutine  RPROP, 
because  in  general  they  will  vary  from  point  to  point  in  the  deposit 
even  though  the  soil  is  homogeneous,  i.e.,  of  one  material  type. 

SIGB:  [olN  j;  i.e.,  the  total  stress  at  the  beginning  of  the  increment; 

compressive  stresses  are  taken  to  be  positive. 

EPB;  [e]^  l5  t^ie  stra*n  at  t^e  beginning  of  the  increment;  compressive 

strains  are  taken  to  be  positive. 

DSlGs  [4a]^^  i.e.,  the  estimate  (supplied  by  the  parent  finite  element 

*  # 
program)  of  the  total  stress  increment. 

DEP:  [Ae]nk  j*  i.e.,  the  estimate  (supplied  by  the  parent  finite  element 

9  * 
program)  of  the  strain  increment. 

DB  &  DE:  lnJN-1  K1  or  [D  1N_1>K_j  and  or  tr>tlN,K-l  (see 

explanation  for  KIND);  i.e^  the  estimates  of  the  incremental  stress- 


A  discussion  of  the  solution  of  initial  estimates  for  these  quantities  is  given 
in  Appendix  D. 


UB: 

DLTAU: 

GAM: 

KIND: 


strain  properties  for  eqs.  (17)  and  (IS)  calculated  by  the  subroutine 
and  supplied  to  the  parent  finite  element  program. 

j;  i.eM  the  pore  water  pressure  at  the  end  of  increment  N-J  (see 
explanation  for  KIND). 

^UN  K  1’  ue’’  est*mate  P°re  water  pressure  increment 

(see  explanation  for  KIND). 

T;  i.e.,  the  combined  bulk  modulus  for  the  soil  particles  and  the  pore 
water  (see  explanation  for  KIND). 

A  flag  assigned  a  value  of  zero  or  one,  depending  on  how  undrained 
conditions  are  being  modeled  in  the  parent  finite  element  program. 
A  value  of  zero  is  required  when  the  special  formulation  for 
incompressible  and  nearly-incompressible  solids  [18,20]  is  being  used 
(i.e.,  the  pore  water  pressure  is  treated  as  a  primary  dependent 
variable  at  the  global  level).  In  such  cases  the  parent  finite  element 
analysis  calculates  u^  ^(UB)  and  Au^  ^  ^  (DLTAU)  at  the  global 
level  and  supplies  them  to  the  subroutine.  The  value  of  F(GAM)  is 
required  at  the  global  level  for  the  near ly-incompressible  formulation, 
and  is  supplied  to  the  parent  program  by  the  subroutine.  The  OB 
array  is  the  [Dj^  j  ^  j  matrix  of  eqs.  (17)  and  (18),  etc. 

A  value  of  one  is  used  for  KIND  when  the  conventional 
formulation  for  compressible  solids  is  used  for  both  drained  and 
undrained  conditions,  [21].  (That  is  the  only  primary  dependent 
variables  are  displacements,  the  pore  water  pressure  is  treated  as  a 
secondary  dependent  variable).  In  such  cases  the  subroutine  calculates 
the  values  of  uKJ  ^  (Ufi)  and  Au^  K  ^  (DLTAH),  and  supplies  them 
to  the  main  program  for  printing  purposes;  the  DB  array  is  the 
[D*!^  j  ^  j  matrix  of  eq.  02),  etc. 
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LARGE:  A  flag  assigned  a  value  of  zero  if  engineering  strains  are  used  and 

a  value  of  one  if  logarithmic  (natural)  strains  are  used, 

LOCIT:  Maximum  number  of  local  iterations  per  call,  recommended  values 

are  5-10. 

TH1:  Value  of  0j  used  in  eq.  (17). 

Subroutine  CLAY  computes  three-dimensional  incremental  properties.  The 
ordering  of  the  stress  and  strain  components  in  the  fa}  and  { e}  vectors  are 
indicated  in  eq.  (10). 

The  subroutine  can  also  be  used  to  supply  properties  for  two-dimensiona: 
finite  element  analyses.  The  procedure  for  its  use  in  such  cases  and  the  value 
of  the  parameter  IDIM  in  the  subroutine  call  are  described  in  the  following 
paragraphs. 

Axisymmetric  Analysis  (IDIM=3):  The  ordering  of  the  stress  and  strain  components 

are  as  follows  (a  ,aQ,a  ,x  Q,  0.0, 0.0)  and  (e  ,eQ,e  ,y  Q, 0.0, 0.0).  The  indicated 
r  9’  7/  r0  r  0  ?/  r0 

zero  values  must  be  supplied  by  the  parent  finite  element  program  in  the  fo}^, 
{e}^,  {Ao}^^  j  and  {Ae}^^  ^  vectors.  The  incremental  properties  of  interest 
are  in  the  upper-left  4x4  corners  of  the  6x6  Dft  and  HE  arrays  returned  by 
the  subroutine. 

Plane  Stress  (1QIM=3):  The  ordering  of  the  stress  and  strain  components  are  as 

follows  (a  ,a  ,0.0,T  ,0.0,0.0)  and  (e  ,e  ,e  ,y  ,0.0, 0.0).  The  subroutine  can  only 

a  j  Ay  a  y  z  xy 

be  used  for  supplying  properties  for  plane  stress  finite  element  analyses  that 
calculate  the  thickness  strain,  ez,  at  the  global  level.  This  limitation  should 
not  be  of  any  consequence,  because  few  earth  structural  problems  are  plane 
stress  in  nature. 

Plane  Strain;  The  subroutine  can  be  used  to  supply  properties  for  plane  strain 
finite  element  analyses  in  two  different  ways.  For  plane  strain  programs  which 
calculate  the  stress  (o2)  normal  to  the  plane  of  the  body,  IDIM  is  given  the 


value  of  3  and  the  ordering  of  the  stress  and  strain  vectors  are 


i 

\ 


(o  ,0  ,o  ,T  ,0.0, 0.0)  and  (e  ,e  ,0.0,y  ,0.0, 0.0),  respectively.  For  plane  strain 
x  y  z  xy  x  y  xy 

analyses  that  do  not  calculate  the  value  of  IH1W  is  given  the  value  of  2 
and  the  stress  and  strain  vectors  are  (o  ,o  ,T  ,0.0,0. 0,0.0)  and 


(€v’g^yv’0*0’0-0»°*0)’  respectively. 

x  y  a  y 

are  appropriately  arranged  in  each 


The  coefficients  *n  the  HB  and  r>F  arrays 
ase. 


2,4  Application  to  homogeneous  tests 

In  the  assessment  of  the  characteristics  of  a  material  mode,  and  in  the 
fitting  of  it  to  experimental  measurements,  a  means  must  be  available  for  using 
the  model  to  predict  the  results  of  simple  homogeneous  tests.  Program  FVAL 
has  been  written  for  this  purpose  [10,171.  EVAL  can  be  used  for  predicting 
the  behavior  of  homogeneous  soil  samples  subjected  to  arbitrary  homogeneous 
stress  and  strain  histories  for  either  drained  or  undrained  conditions. 

The  solution  history  is  broken  into  'tiistory  segments."  Within  each  history 
segment,  a  consistent  combination  of  six  stress  and  strain  components  are 
prescribed  (i.e.,  the  histones  of  £  or  c  e  or  a  ....and  y  or  t  ).  A 

x  xy  y  yz.  y  z 

different  combination  may  be  prescribed  in  each  history  segment.  For  example, 
a  uniaxial  test  might  involve  two  segments  with  a  specified  value  of  axial  strain 
achieved  at  the  end  of  the  first  segment  and  with  unloading  to  zero  axial  stress 
specified  in  the  second.  Each  history  segment  is  broken  into  increments  with 
iteration  conducted  within  each  increment. 

The  analysis  conducted  by  EVAL  is  essentially  a  one-element  finite  element 
analysis  of  a  homogeneous  body.  For  illustrative  purposes  the  user  can  choose 
to  perform  either  a  conventional  or  "reformulated"  analysis  (see  Section  2.3). 
When  both  analyses  converge,  they  give  identical  results.  For  undrained  and 
near-failure  conditions,  the  reformulated  analysis  will,  in  certain  cases,  converge 
when  the  conventional  analysis  will  not. 
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The  reformulated  analysis  is  a  modification  of  the  mixed  finite  element 
procedure  reported  in  [18].  The  strain  components  are  augmented  by  the  pore 
water  pressure  to  form  the  "mixed”  set  of  primary  dependent  variables.  The 
set  of  governing  equations  are  made  up  of  the  incremental  effective  stress-strain 
relations  (eq.  (16))  and  the  expression  0 =T  (Ae ^  +  Aey  +  Ae^)  -  Au. 

Because  of  the  well  behaved  numerical  characteristics  of  the  model  and 
the  simplicity  of  the  analysis  for  homogeneous  tests,  the  method  of  successive 
approximation  was  found  to  be  entirely  adequate,  i.e,,  8^  =  1 /2.  The  analysis  is 
straightforward  and  well  behaved  with  only  two  special  features  worth  noting. 

Because  the  analysis  can  be  used  for  the  extreme  cases  when  either  all 

the  stress  or  all  the  strain  components  are  specified,  both  the  stress  and  strain 

vectors  are  checked  for  convergence.  The  convergence  check  on  the  stress 

increment,  however,  must  be  done  with  some  care.  The  problem  is  that  the 

relative  measure  of  error,  L^Ao^-Acr^  used  for  the  check 

is  meaningless  when  applied  for  near  failure  conditions  (because  Ao.  =  0, 

Lj(Aa^^)=0).  To  avoid  this  difficulty,  the  denominator  of  the  error  measure 

is  limited  to  a  minimum  value  of  L.(a.r  ,)/10.  The  L,  norm  is  the  sum  of 

IN-1  1 

absolute  values. 

The  second  feature  involves  the  starting  estimates  for  the  strain  increment; 
the  procedure  outlined  in  Appendix  n  is  used. 

The  "input"  instructions  for  EVAL  are  given  in  Appendix  I  and  a  program 
listing  in  Appendix  IV. 
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Appendix  1:  Input  Instructions  tor  EVAL 

I.  Heading  information 
Line  i  (OOA2) 

Columns 


1  -  80 

ITITLE  s 

Information  that  is  to  be  printed  as  a  title 
for  the  analysis 

Initial  State  Parameters 

Line  1  (3E10.3) 

Columns 

1  -  10 

%  ! 

Initial  value  of  effective  preconsolidation 
pressure 

11  -  20 

'o  ! 

Initial  void  ratio 

21-30 

%  ! 

# 

Initial  hydrostatic  confining  pressure 

The  strains  produced  by  the  application  of  are  not  calculated.  It  is 
assumed  that  o£  is  applied  under  drained  conditions  (regardless  of  value  of 
D.  If  it  is  desi&d  to  calculate  the  strains  due  to  oc  and/or  to  apply  0C 
under  undrained  conditions,  then  oc  is  set  equal  to  zero  and  the  confining 
pressure  is  applied  in  history  segment  1. 
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UL  Material  Parameters 


Line  1  (8E10.3) 


Columns 


1  -  10  X 

11-20  k 


Slope  of  isotropic  consolidation  line  for  an 
e-tn  p'  plot 

Slope  of  elastic  rebound  line  for  an  e-£n  p' 
plot 


21-30  M 

c 


Slope  of  critical  state  line  in  triaxial  space 
(compression) 


31-40  R 

c 

41-50  A 

c 

51-60  T 


Parameters  describing  shape  of  bounding 
surface  (compression) 


61  -  70 


71-80  v  or  G 


Transitional  value  of  confining  pressure 

separating  linear  rebound  curves  on  e-fcn  p> 

and  e-p'  plots.  Suggested  range  of  values  = 

*3P  +  1.0P 
a  a 

Poisson’s  ratio  or  shear  modulus 


Line  2  (8E1Q.3) 


Columns 


1  -  10 

P 

a 

11  -  20 

r 

21 

-  30 

m 

31 

-  40 

h 

c 

41 

-  50 

h2 

51 

-  60 

n=M  /M 
e  c 

61 

-  70 

“=Vc 

71 

o 

00 

1 

r=Re/Rc 

Atmospheric  pressure 

Combined  bulk  modulus  for  soil  particles  and 
pore  water  (r=0  for  drained  conditions);  if 
no  information  is  available,  it  is  suggested 
that  a  value  of  20,000  P  be  used 

Hardening  parameter 

Shape  hardening  parameter  for  compression 
Shape  hardening  parameter  on  the  I-axis 


Ratio  of  extension  to  compression  values 


Note:  The  input  in  this  group  is  read  by  subroutine  RPROP,  detailed  definitions 
of  the  several  material  parameters  are  given  in  [5,13  ]. 
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Line  3  (3E10.3) 
Columns 


1  -  10  a=Ae/Ac 

11-20  C 

21-30  s 


Ratio  of  extension  to  compression  values 

Projection  center  variable 

Elastic  zone  variable  (a  value  of  1.0  gives  no 
elastic  zone) 


IV.  Iteration  Information 

Line  1  (I5.2E10.2.3I5) 

Columns 

1  -  3  ITMAX  :  Max.  no.  of  iterations  per  increment  (values 

of  5-10  are  suggested) 

6-15  ERMAX  :  Maximum  permitted  relative  difference 

AL^/Ly.  for  the  norms  of  the  incremental 
strAs  and  strain  vectors.  Values  of  .01  to 
.001  are  suggested.  If  convergence  does  not 
occur  in  ITMAX  iterations,  the  program  prints 
a  message  and  then  continues  to  the  next 
problem 

16-25  CONFL  :  (0.0  <  (  )  <  1.0)  Establishes  upper  and  lower 

limits  ofl/CONFL  and  CONFL  for  the 
calculated  values  of  the  Aitken's  acceleration 
factors  (if  it  is  desired  not  to  use  acceleration 
factors  set  CONFL  =  1.0). 


26  -  30  KINO 


0  Reformulated  (for  nearly -incompressible 

conditions)  analysis 

1  Conventional  analysis 


1 0 


31  -  35 


LARGE  : 


1 


Engineering  stresses  and  strains  used  in 
analysis 

True  stresses  and  strains  used  in  analysis 


36-40 


LOCIT 


Maximum  number  of  iterations  to  be  used 
locally  within  subroutine  CLAY  (default=l) 
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V.  Output  Control  Information 


Line  1  (215) 
Columns 

1  -  5 


IPRNTl 


6-10 


IPRNT2 


0  print  incremental  stress  and  strain 
values 

1  suppress  printing  of  incremental 
values 

0  print  computed  incremental  critical 
state  stresses  and  strains 

1  suppress  printing  of  incremental 
critical  state  values 


VI.  Printer  Plot  Control  Information 

Line  1  10(»X,I1) 

Columns 


51 


1PLOT  (I) 


c 


Do  not  generate  printer  plot  type  I 

* 

Generate  printer  plot  type  I 


VO.  Description  of  M  history  segments 


For  each  of  the  M  history  segments,  one  line  (6(I1,E9.3),  15,  El 0.3)  is 
required: 


Columns 


ICOD 


0  -  a 


M 


1  -  € 


2  -  10 


M 


value  of 


is  specified 


M 


for  1C,  = 


-i 


*  For  1=1,6  —  for  an  explanation  of  the  contents  of  plot  ”1"  see  note  following 
"experimental  data  for  plotting"  (section  IX). 

** 

0x  is  the  value  of  a  at  the  end  of  segment  W,  etc.  Thus,  when  IC^  = 

0,  ^  change  in  of  (a  -a  ^  iS  aPP^e<^  in  NINC  increments  during 

the  loading  segment,  a ^  is  value  of  ox  calculated  (ICj  =1)  or 

specified  (IC.  =0)  at  th^'end  of  segment  M-l. 

M-l 
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61  -  65 


NINC 


number  of  increments  into  which  segment 
W  is  to  be  divided 


66  -  75  SR  =  Increment  ratio  for  specified  quantities 

(e.g.,  if  IC! |=0  then  Aax  /Aox 
SR;  A  a  is  the  increment  oY  a  a^pfied 

x\j  x 

during  increment  N  of  loading  segment  M, 
etc.).  A  value  of  1.0  gives  equal 
magnitude  increments  for  the  segment. 

VIIL  Termination  of  Loading  History  Input 

Line  1  (II) 

Column 

1  :  enter  the  integer  9  to  signify  the  end  of 

loading  history  input 

DC.  Experimental  data  for  plotting 

F°r  each  of  the  "printer  plots”  requested  on  the  ’’printer  plot  control  line” 

the  following  information  must  be  specified. 

Line  1  (15) 

Columns 

1-5  NE  =  Number  of  experimental  points  to  be 

printer  plotted  (symbol  //).  The  calculated 
points  will  likewise  be  printer  plotted 
(symbol  *). 

Lines  2-(N  +  l)  (8E10.3),  where  N  .  NE/S 

Columns 

1-10 

el 

11-20  xo 

e2 

21  -  30  xo 

e3 

31  -  40  x 

e4 

41  -  50  x 

e5 

51  -  60  x 

e6 

61  -  70  x 

e7 

71  -  SO  x 

es 


Lines  2  contains  the  abscissas  of  the  first 
eight  experimental  points,  etc. 
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Columns 


1-10  y_ 

el 

11  -  20  y. 

e2 

21-30  y. 

e3 

31  _  40  y  Line  (N+2)  contains  the  ordinates  of  the 

e4  first  eight  experimental  points,  etc. 

41-50  y. 

e5 

51-60  y 

e6 

61  -  70  y_ 

e7 

71  -  80  y 

e8 

Note:  Currently  6  plotting  slots  (the  maximum  number  of  plots  can  easily 

-  ■  —  be  increased)  are  in  use  (the  items  being  plotted  in  a  given  plot  can 

easily  be  changed),  i.e. 

1  =  1:  Q  =  (oz  -  ox)  vs.  P  =  ((dx  +  oy  +  az)l 3  -  u) 

1  =  2:  Q  vs.  ez 
1  =  3:  u  vs.  ez 

I  =  4:  AV/V  (volume  change)  vs.  e 
o  z 

I  =  5:  Q/Po  vs.  P/PQ 
1  =  6:  Q/Pq  vs.  ez 


The  above  input  sequence  (sections  I-IX)  is  repeated  for  each  subsequent  analysis. 


Appendix  II:  Initial  Estimates  for  the  Stress  and  Strain  Increments 


The  matrices  [D 1  ^  j  ^  j  and  |  are  calculated  using  the  K-l 

est'mates  of  the  stress  and  strain  vectors,  thus  at  the  beginning  of  the  iteration 
process,  initial  estimates  are  required.  For  the  first  iteration  of  the  f.rst 
increment,  they  are  usually  taken  to  be  zero.  For  the  intial  iterations  of 
succeeding  increments,  they  also  can  be  started  at  zero;  however,  it  is  usually 
desirable  to  make  use  of  information  from  the  previous  increment  to  obtain 
better  starting  values.  The  simplest  procedure  is  to  use  as  the  intial  estimate 
the  final  values  found  in  the  previous  increment  multiplied  by  the  ratio  of  the 
time  increments  involved.  This  practice  is  based  on  the  assumption  of  relatively 
uniform  behavior  from  increment  to  increment.  Difficulties  can  arise  when  the 
histories  of  applied  external  loads  or  displacements  acting  on  the  structure  cause 
a  switch  from  loading  to  unloading  in  an  unstable  material  response  regime. 
For  example,  consider  the  one-dimensional  response  shown  in  Figure  3.  Consider 
the  case  when  the  state  of  the  soil  is  at  point  "A"  at  the  end  of  increment 
N-l.  If  during  increment  N,  Ao^  is  specified,  two  final  states  B  and  B'  are 
possible.  One  corresponds  to  Ae  (negative)  and  the  other  to  Ae'  (positive). 
Without  any  additional  information,  no  choice  can  be  made  between  B  and  B'. 
(It  is  easily  seen  that  for  specified  stress  increments  in  the  stable  region  of 
behavior  and  for  specified  strain  increments  anywhere,  no  such  problems  exists.) 
The  suggested  solution  to  this  impasse  is  to  assume  that  the  user  would  not 
attempt  a  stress  controlled  specification  for  "loading"  conditions  (path  A-B')  in 
an  unstable  region  and,  hence,  if  the  stress  increment  is  specified,  unloading  is 
the  proper  behavior  (path  A-B)*.  For  stress  controlled  conditions,  the  selection 
of  the  unloading  path  can  be  assured  if  the  starting  estimate  of  strain  is  of 
opposite  sign  to  that  calculated  in  the  previous  increment.  Thus,  the  following 

*  This  argument  requires  that  the  arrival  at  A  must  have  involved  strain 
controlled  steps. 
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STRAIN  € 

Fipn  3.>  One- Dimensional  Soil  Response 
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strategy  is  recommended.  When  considering  a  series  of  increments  for  which 

11 

the  rates  of  the  externally  applied  loads  and  displacements  do  not  change  sign, 
and  { Ao}^  j  are  used  as  starting  estimates  for  increment  N.  However, 
as  one  such  solution  history  segment  is  ended  and  a  new  one  begins,  the 
prerequisite  conditions  for  the  non-uniqueness  problem  may  occur.  Hence,  for 
the  first  increment  of  each  such  series,  it  is  suggested  that  the  starting  strain 
estimate  be  taken  as  some  small  negative  multiple  (e.g.  -.01)  of  the  value  found 
in  the  previous  increment  (the  stress  increment  would  be  used  unchanged).  The 
reduction  in  absolute  magnitude  is  in  deference  to  the  greater  stiffness 
encountered  in  unloading.  Such  an  initial  estimate  will  force  the  solution  to 
select  path  A  -*■  B  if  the  necessary  conditions  exist  for  the  non-uniqueness  to 
occur.  If  non-uniqueness  is  not  a  problem,  the  only  effect  of  this  procedure  is 
to  slightly  slow  the  convergence  process. 


Xi, 

It  is  assumed  that  this  condition  is  sufficient  to  prevent  a  general  switch 
from  loading  to  unloading  within  the  soil  mass. 
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SUBROUTINE  CLAY  dDIM, INC, ITNO, ERMAX, PROP, STOR, SIGBM.EPM, 

•  DSIGM , DEPM ,DB,DE, UB, DLTAU , GAMMA , KIND , LARGE , LOCIT, TH 1 ) 

C  IHIHIIHHIIIHIIIHIHIMHIIIHIIHtlHIlHlllimillHIIMIIII 

C  *  Subroutine  to  evaluate  Iannis  Dafalias'  bounding 

C  1  surface  plasticity  model  for  clay  soils. 

C  *  Fortran  77  version,  Prepared  by  L.R.  Herrmann  and  V.Kaliakln 

C  •  at  the  University  of  California,  Davis  Campus. 

C  eaaaeeeeaeeeeeeeeeeeeeeeeeeaeeeeeeaaaiaeeeeeeeeeeeeeeeeeeeeaaaaeate 

INTEGBR  I, J,K, IT, IDIH, INC, ITNO, KIND, LARGE, LOCIT, 11(6) 

REAL  PROP ( 1 9 ) , STOR ( 7 ) , SIGBM ( 6 ) , EPM ( 6 ) , DSIGM ( 6 ) , DEPM ( 6 ) , DB ( 6 , 6 ) , 

•  DB(6,6),SIGB(6),EPB(6),DSIO(6)1DEP(6),DEPT(3,3),SB(3,3), 

•  SB(3, 3), DLTA(3, 3), ERMAX, UB, DLTAU, GAMMA, GAMMAT, SMALL, DIL, 

•  DDIL , VOIDB , VOIDE , IIB , HE , X JB , X JE , SCUBSB , SCUBEE , SIN3AB , 

•  SIN3 AE , COS3 AB , C0S3AE , BULKB , BULKE , GB , GE , XIOB , IIOE , XIL , 

•  TH 1 , TEMPI , TEMP2 , TEMP3 , TEMP* 

C 

DATA  11/11,22,33,12, 13,23/,  DLTA/1 .0,3*0. 0, 1 .0,3*0. 0, 1 .0/ 
SMALL=0.0001*PR0P(8) 

C 

DO  100  1=1,6 

SIGB( I )=SIGBM( I ) 

DSIG ( I ) =DSIGM ( I ) 

EPB(I)=EPM(I) 

DEP(I)=DEPM(I) 

100  CONTINUE 

Initialize  history  if  necessary 

IF(INC  .EQ.  1  .AND.  ITNO  .EQ.  1)  THEN 
STOR ( 2 ) =ST0R ( 1 ) 

ST0R(3)=0.0 

STOR(5)=0.5*(SIGB(1)  SIGB(2)) 

STOR(6)=0.01*PROP(8) 

ELSE 


Update  history  if  necessary 

IF(INC  .GT.  1  .AND.  ITNO  .EQ.  H  THEN 
STOR ( 1 ) =STOR ( 2 ) 

STOR(3)=STOR(3)  +  ST0R(4) 

STOR ( 5 ) =ST0R ( 5 )  +  ST0R(6) 

END  IF 
END  IF 


Convert  from  plane  strain  to  3-dimensonal  state  if  necessary 
IF(IDIM  .EQ.  2) 

•  CALL  TVODIM  (1 ,SIGB,EPB, DSIG, DEP,DB,DE, STOR) 


DIL=0.0 
DDILsO.O 
DO  200  1=1,3 

DIL  =DIL  ♦  EPB(I) 

DDILsDDIL  ♦  DEP(I) 

200  CONTINUE 

Determine  3-dimensional  incremental  properties. 
Iterate  on  the  stress  estimate. 

DO  600  IT«1, LOCIT 


38 


ooo  ooooo  ooo  ooo  ooo 


Compute  values  of  the  Invariants 
IPCIT  .EQ.  1) 

•  CALL  INVAR  ( 1 ,  PROP, STOR,SIGB,DSIG, DEP, DEPT, VOIDB, KIND, 

•  LARGE , SMALL , XIB , XJB , DIL , DDIL , SB , SCUBEB, SIN3 AB , C0S3AB , 

•  GAMMA, GAMMAT,UB,DLTAU, II, DLTA) 

CALL  INVAR  (2, PROP, STOR,SIGB,DSIG, DEP, DEPT, VOIDE, KIND, 

•  LARGE , SMALL , XIE , X JE , DIL , DDIL , SE , SCUBEE , SIN3AE , COS3 AE , 

•  GAMMA , GAMMAT , UB , DLTAU , II , DLTA ) 

Calculate  elastic  incremental  properties 
IF(IT  .EQ.  1) 

•  CALL  ELASTC  ( PROP, VOIDB, XIB, DB, BULKB, GB, GAMMAT, II , DLTA ) 
CALL  ELASTC  (PROP, VOIDE, XIE, DE.BULKE.GE, GAMMAT, II, DLTA) 

Calculate  the  size  of  bounding  surface 

XIOB=STOR( 1 ) 

XI0E=ST0R(2) 

XIL=PR0P(7 ) 

TEMPI = 1. 0/ ( PROP( 1 )  -  PR0P(2) ) 

TEMP2= (XIE  -  XIB)/3.0 

IF(XIOB  .GE.  XIL  .AND.  XIOE  .GE.  XIL)  THEN 

XI0E=XI0B,EXP(TEMP1*0.5#( (VOIDB  +  VOIDE)»DDIL 

•  -  (VOIDB/BULKB  +  V0IDE/BULKE)»TEMP2) ) 

ELSE 

TEMP3=XIOB 

IF(XIOB  .LT.  XIL)  TEMP3=XIL 
TEMP4=XIOE 

IF(XIOE  .LT.  XIL)  TEMP4=XIL 

XIOE=XIOB  +  TEMP 1*0. 5* ( (TEMP4*VOIDE+TEMP3#VOIDB)*DDIL 

•  -  (TEMP3*V01DB/BULKB  +  TEMP4«VOIDE/BULKE)«TEMP2) 

END  IF 

ST0R(2)=XI0E 

IF(INC  +  ITNO  +  IT  .GT.  3)  THEN 


Calculate  the  bounding  surface  parameters  and  parameters  associated 
with  the  plastic  portion  of  the  incremental  properties. 

IF( IT  .EQ.  1  .OR.  (INC  +  ITNO  +  IT  .EQ.  4)) 

•  CALL  PLASTC  (PROP, DEPT, VOIDB, XIB, XJB, XIOB, DDIL, 

•  SB, SCUBEB, SIN3AB, COS3AB, DB, BULKB, GB, II , DLTA ) 
CALL  PLASTC  (PROP, DEPT, VOIDE, XIE, XJE,XIOE, DDIL, SE, SCUBEE, 

•  SIN3AE.COS3AE, DE.BULKE.GE, II, DLTA) 

Calculate  revised  total  stress  estimates  and  error  norms 
IF(IDIM  .EQ.  2) 

•  CALL  TWODIM  (2,SIGB,EPB,DSIG,DEP,DB,DE,ST0R) 

TEMP2=0. 0 

TEMP3=0.0 
TEMP4=0. 0 
DO  400  1=1,6 
J=II(I)/10 
K=MOD(II(I ) , 10 ) 

TEMP 1 =DLTA (K, J) *DLTAU 
DO  300  J=1 ,6 
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TEMPI sTEMPI  ♦  ((1.0  -  TH1 )*DB(I, J) 

•  ♦  TH1*D£(I, J) )*DEP(J) 

300  CONTINUE 

TBMP2*TEMP2  ♦  ABS( TEMPI  -  DSIG(I)) 

TEMP3*TEMP3  ♦  ABS( TEMPI) 

TEMP4s TEMP4  ♦  ABS(SIGB(I)) 

DSIG(I). TEMPI 
400  CONTINUE 

IF(TEMP3  -LT.  TEMP4«0.01)  TEMP3«0. 01 »TEMP4 
IF(TEMP3  .HE.  0.0)  THEN 

IF(TEMP2/TEMP3  .LT.  ERMAX)  GOTO  700 
END  IF 
END  IF 

600  CONTINUE 
700  CONTINUE 

Convert  3-dieensional  properties  to  plane  strain  if  necessary 
IF(IDIM  .EQ.  2) 

•  CALL  TWODIM  (3,SIGB,EPB,DSIG,DBP,DB,DE,STOR) 

RETURN 

END 
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SUBROUTINE  INVAR  (IK, PROP, STOR,SIG,DSIG, DEP, DEPT, VOID, KIND, 

•  LARGE , SHALL , XI , X J , DIL , DDIL , S , SCUBE , SIN3 A , C0S3 A , 

•  GAMMA, GAMMAT,UB,DLTAU, II, DLTA) 


Subroutine  to  compute  values  of  Invariants 
INTEGER  I, J,K,N, IK, KIND, LARGE, 11(6) 

REAL  PROP(19),STOR(7),SIG(6),DSIG(6),DEP(6),DEPT(3,3),S(3,3), 

•  DLTA (3,3), VOID , SMALL , XI , X J , DIL , DDIL , SCUBE , SIN3A , C0S3A , GAMMA , 

•  GAMMAT , UB , DLTAU , ARB , FACTOR , TEMP 1 

DATA  ARB/ 1000.0/ 

FACTORrO.O 

IF (IK  .GT.  1)  FACT0R= 1 . 0 

XI  =0.0 

XJ  =0.0 

SCUBE=0.0 

SIN3A=0.0 

Calculate  first  stress  invariant 
DO  100  1=1,3 

XI=XI  +  SIG(I)  +  FACTOR*DSIG(I ) 

100  CONTINUE 

V0ID= 1.0  +  ST0R(7) 

IF(LARGE  .NE.  0) 

•  VOID=VOID*EXP(-DIL  -  FACTOR*DDIL) 

Change  tensor  components  to  matrix  components, 
calculate  deviatoric  stresses. 

DO  200  N= 1 ,6 
I=II(N)/10 
J=M0D(II(N) ,10) 

S(I, J)=SIG(N)  +  FACT0R*DSIG ( N )  -  XI»DLTA (I , J)/3- 0 
S(J,I)=S(I,J) 

DEPT ( I , J ) =DEP ( N ) • ( 1 . 0  +  DLTA(I, J) )*0.5 
DEPT ( J , I ) =DEPT ( I , J  ) 

200  CONTINUE 


Convert  total  stresses  to  effective  stresses 

GAMMA=PR0P(6) 

GAMMAT=0. 0 
IF(KIND  .NE.  0)  THEN 
GAMMAT=GAMMA 
UB=ST0R(3) 

DLTAU=GAMMA»DDIL 
END  IF 

XI=XI  -  3.0»(UB  +  FACTOR»DLTAU ) 

ST0R(4 )=DLTAU 


Avoid  near  zero  value  of  the  first  stress  invariant 

IF(ABS(XI)  .LE.  SMALL)  THEN 
TEMPI =XI 
XIsSMALL 

IF(TEMP1  .LT.  0.0)  XI=-SMALL 
END  IF 
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Compute  the  square  root  of  the  second  deviatoric  stress  invariant 
as  well  as  the  third  deviatoric  stress  invariant 

00  300  1*1,3 

DO  300  J*1 ,3 

XJ=XJ  ♦  S(I,J)«S(I,J) 

DO  300  K=1,3 

SCUBEsSCUBE  +  S(I, J)*S(J,K)«S(K,I) 

300  CONTINUE 

SCUBE=SCUBE/3 . 0 

Arbitrary  check  to  avoid  excessively  snail  values  of  J 

X J=SQRT ( 0 . 5#X J ) 

IF(XJ*ARB  .LT.  XI)  XJ=0.0 

Compute  the  sine  and  cosine  of  three  times  the  "Lode"  angle 

IF(XJ  .GT.  SMALL)  SIN3A=1 .5«SQRT(3.0)«SCUBE/XJ”3 
IF(SIN3A  .GT.  1.0)  SIN3A=  1.0 
IF(SIN3A  .LT.  -1.0)  SIN3A=-1.0 

COS3A=SQRT(1.0  -  SIN3A»»2) 

RETURN 

END 
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SUBROUTINE  ELASTC  (PROP, VOID, XI, D, BULK, G.GAMMAT, II, DLTA) 

C 

INTEGER  I,J,K,L,M,N,II(6) 

REAL  PROP(19),D(6,6),DLTA(3,3),VOID,XI,BULK,G,GAMMAT,TEMP1, 

•  TEMP2,TEMP3 

Calculate  the  bulk  and  shear  moduli 

TEMPIsl .5*0 .0  -  2.0»PROP(5))/(1.0  +  PR0P(5)) 

TEMP2= VOID/3 . O/PROP ( 2 ) 

TEMP3=XI 

IF(TEMP3  .LT.  PR0P(7 ) )  TEMP3=PROP(7 ) 

BULK=TEMP2»TEMP3 
IF(PR0P(5)  .LE.  0.5)  THEN 
G=TEHP1  •  BULK 

ELSE 

G=PR0P(5 ) 

END  IF 

Calculate  elastic  incremental  properties 

DO  100  M=1,6 
I=II (M)/10 
J=MOD(II(M), 10) 

DO  100  N=M,6 
K=II(N)/10 
L=M0D(II(N) ,10) 

TEMP 1 =DLTA ( K , I ) *DLTA ( L , J )  +  DLTA ( K , J ) *DLTA ( I , L ) 
D(M,N)=TEMP1*G  +  (BULK  +  GAMMAT 

•  -  2 . 0*G/3 . 0 ) *DLTA ( I , J ) *DLTA ( K , L ) 
D(N,M)=D(M,N) 

100  CONTINUE 
RETURN 
END 
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SUBROUTINE  PLASTC  (PROP, DEPT, VOID, XI, XJ, XIO, DDIL, S, SCUBE, 

•  SIN3A , C0S3A , D , BULK , G , II , DLTA ) 

INTEGER  I,J,K,L,M,N,LL,LFLAG,II(6) 

REAL  PR0P(19),DEPT(3,3),S(3,3),D(6,6),DLTA(3,3),V0ID,X,XI, 

•  X J , XIO , DDIL , SCUBE , SIN3A , C0S3A , BULK , G , XKP , XKPB , BETA , 

•  DBETA , DPI , DF J , DFAL , DF J J , TEMP , TEMP 1 , TEMP2 , TEMP3 , TBMP4 , 

•  TEHP5 , TEMP6 

Calculate  bounding  surface  parameters 

CALL  BOUND  (PROP, VOID,!, XI, XJ, XIO, SIN3A, XKPB, BETA, 

•  DFI , DFJ , DFAL , DFJ J ) 

DBETAsBETA  -  1.0 

IF (DBETA  .LT.  0.0)  DBETA=0. 0 


Check  for  elastic  zone 


LFLAGsO 

TEMPI sBETA  -  DBETA»PROP( 15 ) 
IF(TEMP1  .GT.  0.0)  THEN 
LFLAG=1 


Calculate  the  plastic  nodulus  and  loading  function 

CALL  LODFUN  ( PROP, DEPT, VOID, X,XJ, XIO, DDIL, S, SCUBE, SIN3A, 

1  COS3A, BULK, G, XKP, XKPB.DBETA, TEMPI , DFI, DFJ, DFAL, DFJJ.LFLAG) 

END  IF 


Calculate  plastic  portion  of  the  increaental  properties 

IF(LPLAG  .NE.  0)  THEN 
DO  200  M=1 ,6 
I=II(M)/10 
J=MOD(II(M), 10) 

DO  200  N=M,6 
K=II(N)/10 
L=MOD(II(N), 10) 

TEMP=0.0 
TEMPI =0.0 
TEMP2=0.0 

TEMP3=3.0»BULK*DFI 
TEMP4=G»DFJJ 
TEMP5=SQRT( 3. 0 )«G»DFAL 

TEMP6=XKP  «■  9. 0*BULK*DFI*DFI  ♦  G»DFJ»DFJ 

•  ♦  G»(DFAL»COS3A)»(DFAL»COS3A) 

IF(XJ«4  .NE.  0.0)  THEN 

DO  100  LL=1 ,3 

TEMPI = TEMPI  ♦  S(I,LL)*S(LL, J) 

TBMP2=TEMP2  «■  S(K,LL)»S(LL,L) 

100  CONTINUE 

TEMPI =TEMP5* (TEMPI /XJ##2  -  1 ,5»SCUBE»S(I,  J)/XJ«4 

•  -2 . 0*DLTA ( I , J ) /3 . 0 ) 

TEMP2=TEMP5* ( TEMP2/X J**2  -  1 .5#SCUBE«S(K,L)/XJ”4 

•  -2.0»DLTA(K,L)/3.0) 

END  IF 

TEMP* ( TEMP3*DLTA ( I , J )  +  TEMP4»S(I,J)  «■  TEMPI) 

•  •(TEMP3*DLTA(K,L)  ♦  TEMP4»S(K,L)  ♦  TEMP2)/TEMP6 

D(M,N)=D(M,N)  -  TEMP 
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SUBROUTINE  LODFUN  ( PROP , DEPT , VOID , X , X J , XIO , DDIL , S , SCOBE , SIN3A , 

1  C0S3 A , BULK , G , XKP , XKPB , DBETA , DDEN , DFI , DF J , DFAL , DFJJ , LF ) 


Subroutine  to  calculate  the  plastic  modulus  and  loading  function 

INTEGER  I , J , K , LF 
REAL  ALFUN,CV,RT,SINV 

REAL  PROP(19),DEPT(3,3),S(3,3),VOID,X,XJ,XIO,DDIL,SCUBE,SIN3A, 

•  C0S3A , BULK , G , XKP , XKPB , DBETA , DDEN, DFI, DF J , DFAL , DF J J , XN , R , 

•  Z , XIL , XM , H 1 , B2 , H , SUM , XLF , TEMP 1 , TEMP2 , TEMP3 , TEMP4 , TEMP5 
ALFUN(CV,RT,SINV)=2.0*RT*CV/(1 .0  +  RT  -  (1.0  -  RT)»SINV) 
XN=ALFUN(PR0P(3),  PR0P(4),  SIN3A) 

R  =ALFUN(PROP(9) ,  PROP( 12) , SIN3A) 

HI = ALFUN ( PROP ( 1 6 ) , PROP ( 1 8 ) , SIN3A ) 

Z=XJ»R/(XN«XIO) 

XIL=PROP(7) 

XM=PR0P(17) 

H2=PR0P( 19 ) 


Calculate  the  plastic  modulus 

(using  modified  formulation  for  continuity  across  the  I-axis) 

TEMP1=1.0/(PR0P(1)  -  PROP(2) ) 

TEMP2=Z«XH 

TEMP3=9 • 0*DFI*DFI  +  DFJ*DFJ/3.0 
TEMP4rXI0 

IF(XI0  .LT.  XIL)  TEMP4=XIL 
H=H1*TEMP2  +  H 2»(1.0  -  TEMP2) 

XKPsXKPB  ♦  H»DBETA/DDEN»TEMP1»VOID«TEMP3»TEMP4 

Calculate  the  loading  function 

SUM=0.0 

TEMP1=0.0 

TEMP2=3.0«BULK«DFI 

TEMP3=G*DFJJ 

TEMP4 =SQRT (3.0) *G*DFAL 

TEMP5=XKP  9.0»BULK»DFI»DFI  +  G*DFJ*DFJ 

•  +  G«(DFAL»COS3A)«(DFAL»COS3A) 

IF(XJ»»2  .NE.  0.0)  THEN 

DO  200  1=1,3 

DO  200  J=1 ,3 

TEMP2=0.0 
DO  100  K=1,3 

TEMP2=TEMP2  ♦  S(I,K)*S(K, J) 

100  CONTINUE 

TEMP1=TEMP1  +  (TEMP2  -  1 . 5»SCUBE»S(I, J)/(XJ«XJ ) ) 

•  •D£PT(I,J)/(XJ»XJ) 

SUM=SUM  +  S(I, J)*DEPT(I, J) 

200  CONTINUE 

TEMPI = TEMPI  -  2.0»DDIL/3.0 
END  IF 

XLF=(TEMP2*DDIL  ♦  TEMP3«SUM  TEMP4«TEMP1 )/TEMP5 
Check  for  unloading 


IF (XLF  .LB.  0.0)  LF=0 

RETURN 

END 
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SUBROUTINE  TVODIM  ( IK , SIGB , EPB , DSIG , DEP , DB , DE , STOR ) 


Subroutine  to  perform  manipulation  of  storage  locations 
for  the  case  of  plane  strain. 


INTEGER  I,  IK 

REAL  SIGB(6 ) , EPB(6 )  ,DSIG(6 ) , DEP(6 ) ,DB(6 , 6 ) ,DE(6 , 6 ) , STOR (7 ) .TEMPI 
C 

IF(IK  .EQ.  1)  THEN 
SIGB(4)=SIGB(3) 

SIGB(3)=STOR(5) 

DSIG(4 )=DSIG(3 ) 

DSIG(3)=STOR(6) 

EPB(4)=EPB(3) 

EPB(3)=0.0 

DEP(4)=DEP(3) 

DEP(3 )=0. 0 
DO  100  1=5,6 

SIGB(I)=0.0 
DSIG(I)=0.0 
EPB(I)  =0.0 
DEP(I)  =0.0 
100  CONTINUE 

ELSE  IF (IK  .EQ.  2)  THEN 
C 

C  Compute  and  store  the  stress  increment  in  the  2-direction 

C 

TEMP1=0.0 
DO  200  1=1,4 

TEMPI =TEMP1  +  0 . 5 • ( DB (3,1 )  +  DE(3 , I ) ) *DEP(I ) 

200  CONTINUE 

STOR(6)=TEMP1 

ELSE 

C 

C  Convert  the  3-dimensional  state  to  one  of  plane  strain 

C 

DO  300  1=1,4 

DB( 3 , I )=DB(4 , I ) 

DB(4 , 1  )  =  0 . 0 
DE(3,I)=DE(4,I) 

DE(4 , I ) =0. 0 
300  CONTINUE 

DO  400  1=1,3 

DB(I, 3 )=DB( 1,4) 

DB(I,4)=0.0 
DE(I , 3 )=DE( I ,4 ) 

DE(I,4 )=0.0 
400  CONTINUE 

END  IF 
RETURN 
END 
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SUBROUTINE  BOUND  ( PROP, VOID, X, XI ,XJ,XIO,SIN3A,XKPB, BETA, 

1  DFI , DFJ , DFAL , DFJ J ) 

Subroutine  to  evaluate  relationship  of  current  stress  state 
to  the  bounding  surface 


INTEGER  IZONE 

REAL  ALFUN, CV,RT,SINV 

REAL  DFUN,FUN,FUNC 

REAL  PROP (19), VOID , X , XI , X J , XIO , SIN3A , XKPB , BETA , DFI , DFJ , DFJ J , DFAL , 

•  XN , DNAL , R , DRAL , A , DAAL , Y , C , ARB , BIG , SMALL , T , Q , QC , QO , FOP , X JO , 

•  BT , RHO , XIBAR , THETA , PSI , GAM , DYSAL , DFOPAL , D JOAL , DBTAL , DRHOAL , 

•  TEMP , TEMP 1 , TEMP2 , TEMP3 , TEMP4 , TEMP5 , TEMP6 , TEMP7 , TEMP8 

DATA  ARB/0.001/,  BIG/1 . 0E+20/,  SMALL/ 1 . 0E-20/ 
ALFUN(CV,RT,SINV)=2.0*RT»CV/(1.0  +  RT  -  (1.0  -  RT)»SINV) 
DFUN(FUN,RT,FUNC)=FUN»»2«(1.0  -  RT)/(2.0»RT»FUNC) 

XN: ALFUN ( PROP ( 3 ) , PROP ( 4 ) , SIN3A ) 

DNAL=DFUN ( XN , PROP ( 4 ) , PROP ( 3 ) ) 

R= ALFUN ( PROP ( 9 ) , PROP ( 1 2 ) , SIN3 A ) 

DRAL=DFUN ( R , PROP ( 1 2 ) , PROP ( 9 ) ) 

A=ALFUN ( PROP (10), PROP ( 1 3 ) , SIN3A ) 

DAAL=DFUN ( A , PROP ( 1 3 ) , PROP ( 1 0 ) ) 

Y=R»A/XN 
C=PROP( 14 ) 


Shift  projection  point 
TEMPI =XI  -  XIO»C 

IF(ABS(TEMP1 )  .LT.  ARB)  TEMP 1= ARB 
TEMP2=C  -  1.0/R 
TEMP3=TEMP1«TEMP2 
TEMP4=C»(C  -  2.0/R) 

Q  =XJ/TEMP1 
QC:XN/( 1.0  -  R»C) 

QO=BIG 

IF(C  .NE.  0.0)  QO=XN»( SORT (1.0  +  Y»Y)  -  (1.0  +  Y))/R/C 
IF(XJ  .EQ.  0.0)  THEN 

IF (TEMPI  .GT.  0.0)  THEN 
IZONE= 1 

ELSE 

IZONE=3 
END  IF 

ELSE  IF(C  .GE.  1.0/R)  THEN 

IF(Q  .GE.  0.0  .OR.  Q  .LE.  QC)  THEN 
IZONE: 1 

ELSE  IF(Q  .GE.  QO)  THEN 
IZONE: 3 

ELSE 

IZ0NE:2 
END  IF 

ELSE 

IF(Q  .GE.  QC)  THEN 
IZONE:2 

ELSE  IF(Q  .GE.  0.0)  THEN 
IZONE: 1 

ELSB  IF(Q  .LB.  QO)  THEN 
IZONE:2 

ELSE  48 
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IZ0NE=3 
END  IF 
END  IF 

IF(IZONE  .EQ.  1)  THEN 
C 

C  Projection  on  ellipse  1 

C 

TEMP5  =  TEMPI  •TEMPI  +  ((If  -  1 . 0)*XJ/XK -**Z 

BETA =X  10* (-TEMP3+SQRT( TEMP 3*TEMP3- TEMPI-. •  ( TEMP4+( 2.  C-R)/R) ) ) 

•  /TEMP5 

ELSE  1F(IZ0NE  .EQ.  2)  THEN 

Projection  on  hyperoo  u 

TEMP5=TEMPt  -  2.0«A, E/XN 
TEMP6=XJ*( 1.0/R  +  A/XN'/XN 
TEMP7=TEMP3  +  TEMP6 
TEMP8=TEMF1 “TEMPI  -  (XJ/ XN)*(XJ/XN> 

BETA=-0.5*XIO»TEMP5/TEMP? 

IF(TEMP8  .NE.  0.0) 

BETAs XI0*(-TEMP7  +  SORT ( TEMP7 “TEMP?  -  TEMP8*TEMP5 ) ) /TEMP6 

ELSE 

Projector!  on  ellipse  2 
?=PR0?( 11) 

TEMP5  =  SQRT( 1 . P  *  Y*Y  ) 

F0P=XN/TEMP5 

XJ0=A*( 1 . 0  +  Y  -  TEMPER 
BT=T*  (XJO  -  T*FOP  .  /  ;  KJO  -  ?.'>.* EOF  ) 

RH0=  ( BT  -  T1/F0P/X..0 
TEMP6=T  -  BT  *  C 

temp7=tempi«:emp  f. 

TEMP8sTEMP1  “TEMPI  -  RH0*X.;*X« 

BETA  =  X 1 0*  f  -TEM i  ?E!  i  ’ 

-TEMP8*  1  Tr.MPn*'.  KM:  •  -  BT*  r  '  .  > 

END  IF 

Compute  derivatives  vif  the  ..fundinp  ;  ;rfaee  w.r.t.  invariants 
and  the  value  of  l.i  r bounding"  p... modulus  for  the 

app--'-...  w-te  zone. 

XIBARsBETA* ( XI  -  XI0»C)  +  XldC 
IF (XIBAR  .EQ.  0.0)  X 7  h A H - S M A 

THETA=BETA*XJ/XLBAit 
X=THETA/XN 
TEMP5-XI0 

IF(XI0  . 1.T.  PROP (7  )  )  TEMPI  =  . 

TEMP=  12.0*V0ID/(  PROP  1  '  •  PhQPd)  ‘.*X  i'r'tX7 f*T  EMP5 

IF (I ZONE  .EQ.  1)  THEN 

Normal  consolidation  .•.■re  I  .  : ■  * 

PSI=Y/((R  -  7.0)*FR  -  1.0)) 

TEMP5=R*( 1 . 0  +  X*X  *  R«(R  -  .0  »X*X‘ 

GAM-  (1.0  +  '  R  -  1.0  »SCJf.M.O  •>  RVR  -  .  . 0  )*X*X  ) )  /TEMP5 
DFI  =  ?.  0#Xi.  -  'GAM  -  .o 

DFJ.J=2.0«XxO*CAM*».  (k  -  •  .0)  XN,  **x*PSi*3ETA/XIBAR 
DFJ=DFJJaXJ 


XKPB=  TEMP* (GAM  -  1.0/R)*(GAM  *  R  -  2. 0)"PSI«PSI/R 
DFAL=PSI"6 . 0" ( R- 1 . 0 ) "THETA "GAM^XIO" ( ( ( R- 1 . 0 ) / ( R"R" 

(2.0/R-GAM-1 .0) )  ♦  1.0)*DRAL  -  (R-1 .0)"DNAL/XN)/(XN"XN) 
RETURN 

ELSE  IF(IZONE  .EQ.  2)  THEN 

Overoonsolidated  zone  (hyperbola) 

TEMP5=1.0  -  X*(1.0  Y) 

GAM=-(TEMP5«-SQRT((X-Y-1.0)""2  +  (X"X-1 . 0)"Y"Y ) )/ (R"(X"X-1 .0) ) 
DFI=2.0"XIO"(GAM  -  1.0/R) 

DFJ=2. 0"XI0"( (1.0  ♦  Y)/R  -  X"GAM)/XN 
DFJJ=DFJ/IJ 

XKPB=TEMP* ( GAM  -  1 .0/R)"(TEMP5"GAM  +  2.0"A/XN)/R 
DFAL=6. 0*110* (DNAL* ( THETA*GAM/XN-1 . 0/R+A/ (R"THETA*GAM) 

-2 . 0"A/XN ) / (XN"XN )+DRAL" ( 1 . 0/THETA- 1 . 0/XN+A/ (XN"THETA"GAM ) ) 
/(R"R)  +  DAAL"( 1 . 0/XN  -  1 . 0/ ( THETA"GAM"R ) ) /XN ) 

RETURN 

ELSE 

Tension  zone  (ellipse  2) 

PSIr1.0/(R"(BT  -  T)) 

GAM= ( -T+BT-SQRT ( BT"BT-RHO"THETA"THETA"T" ( T-2 . 0*BT ) ) ) 

/(1.0  +  RHO"THETA"THETA ) 

DFI=2.0"PSI"XI0"(GAM  +  T  -  BT) 

DFJJ=2. 0*PSI*XIO"GAM*RHO"BETA/XIBAR 
DFJ=DFJJ"XJ 

XKPB-TEMP*PSI*PSI * ( GAM+T-BT )*( CAM* ( BT-T )  +  T"(2. 0"BT-T) ) 
DYSAL=Y"(DRAL/R  +  DAAL/A  -  DNAL/XN) 

DFOPAL=FOP" ( DNAL/XN  -  Y"DYSAL/(1.0  +  Y"Y ) ) 

DJOAL=XJO" (DAAL/A  -  DYSAL/Y)  A"(1.0/Y  -  FOP/XN)"DYSAL 
DBTAL=((T-BT)"DJOAL  -  ( T-2. 0"BT)"T"DF0PAL)/(X JO-2. 0"T"FOP) 
DRHOAL=DBTAL/FOP/XJO  -  RH0"(DF0PAL/F0P  +  DJOAL/XJO) 
DFAL=3.0"PSI"XIO"THETA"GAM"(DRH0AL  +■  2.0"RH0*DBTAL 
/(T>GAM-2.0"BT)) 

RETURN 
END  IF 
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SUBROUTINE  JPLOT  (NP,PINC,IX,IY,XT, YT) 


Printer  plotting  subroutine  written  by  J.  S.  De  Natale. 
Modified  by  V.  Kaliakin,  1983. 

CHARACTERS  XT(8) fYT(8) ,BUF( 128) 

INTEGER  ICPX( 150 ) , ICPY( 150) , 1CEX(20), ICEY(20) , I , J, IX , IY , IK , 

•  JK,NP,NEXP 

REAL  PINC ( 1 50 , 8 ) , XEXP ( 20 ) , YEXP ( 20 ) , XLAB (11), YLAB ( 1 1 ) , XMIN , XMAX , 

•  YMIN, YMAX ,XDIST, YDIST,XINC , YINC , TEMP 

Read  values  for  the  experimental  data  points 

READ (5, 902)  NEXP 
IF ( NEXP  .NE.  0)  THEN 

READ (5, 904  )  (XEXP(I), 1=1, NEXP) 

READ (5, 904 )  (YEXP(I ), 1=1 , NEXP) 

END  IF 


Establish  minimum  and  maximum  axes  values 


XMIN=PINC ( 1 , IX ) 

XMAX=PINC( 1 , IX) 

YMIN=PINC( 1 ,IY) 

YMAX=PINC( 1 , IY) 

DO  100  1=2, NP 

IF (XMIN  .GT.  PINC(I , IX) ) 
IF (XMAX  .LT.  PINC (I , IX ) ) 
IF(YMIN  .GT.  PINC (I , IY ) ) 
IF ( YMAX  .LT.  PINC(I.IY)) 
100  CONTINUE 

IF (NEXP  .NE.  0)  THEN 
DO  200  1=1, NEXP 

IF (XMIN  .GT.  XEXP(I)) 
IF (XMAX  .LT.  XEXP(I)) 
IF(YMIN  .GT.  YEXP(I)) 
IF (YMAX  .LT.  YEXP(I) ) 
200  CONTINUE 
END  IF 


XMIN=PINC(I,IX) 
XMAX=PINC(I,IX) 
YMIN=PINC ( I , IY ) 
YMAX=PINC(I, IY ) 


XMIN=XEXP(I) 

XMAX=XEXP(I) 

YMIN=YEXP(I) 

YMAX=YEXP(I) 


Establish  adjusted  maximum  and  minimum  axes  values 

CALL  SHIFT  (XMIN, XMAX) 

CALL  SHIFT  (YMIN, YMAX) 

Establish  axes  labels 

XDIST=XMAX-XMIN 
YDIST=YMAX-YMIN 
XINC=XDIST/ 10.0 
YINC=YDIST/10.0 
XLAB(  1 )=XMIN 
XLAB(11)=XMAX 
YLAB(  i )=YMIN 
YLAB(11 )=YMAX 

DO  300  1=2,10 

TEMP=FLOAT ( I- 1 ) 

XLAB(I)=XMIN  +  XINC»TEMP 
YLAB(I )=YMIN  YINC»TEMP 
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300  CONTINUE 


Establish  data  point  coordinates 

DO  400  1=1, NP  I  For  computed  data  points 

TBMPs ( PINC ( I , IX )  -  XMIN)/XDIST« 100.0 
CALL  SETCRD  (I.TEMP.ICPX) 

TEHP=ABS( PINC ( I , IT )  -  IMAX)/IDIST«  50.0 
CALL  SETCRD  (I.TEMP.ICPI) 

400  CONTINUE 

IF(NEXP  .NE.  0)  THEN  I  For  experimental  data  points 

DO  500  I=1,NEXP 

TEMPs (XEXP(I )  -  XMIN)/XDIST« 100.0 
CALL  SETCRD  ( I , TEMP , ICEI ) 

TEMP=ABS(IEXP(I)  -  IMAX)/IDIST»50.0 
CALL  SETCRD  (I , TEMP, ICEI) 

500  CONTINUE 
END  IF 

Printer  plot  the  data  points 

JK=0 

WRITE (6, 900) 

DO  800  IK=0,50  I  Loop  through  each  of  the  print  lines 

CALL  BZER0  (BUF.128) 

DO  600  J=1,NP  I  Print  computed  data  points 

IF(ICPKJ)  .BQ.  IK)  THEN 
BUF(ICPX(J)  ♦  25)='*' 

END  IF 

600  CONTINUE 

IF(NEXP  .NE.  0)  THEN  I  Print  experimental  data  points 

DO  700  J= 1 , NEXP 

IF(ICEKJ)  .EQ.  IK)  THEN 
BUF(ICEX(J)  ♦  25)='#' 

END  IF 

700  CONTINUE 

END  IF 

CALL  BORDER  (JK.IK.XLAB.ILAB.XT.IT.BUF) 

800  CONTINUE 

900  F0RMAT( 1H1 ,/) 

902  FORMAT (15) 

904  F0RMAT(8E10.3) 

RETURN 

END 
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SUBROUTINE  SHIFT  (PMIN.PMAX) 


Subroutine  to  compute  adjusted  maximum  amd  minimum  axes  values. 
Written  by  J.  S.  De  Natale. 


INTEGER  II, JJ,NN 

REAL  PMIN,PMAX,DX,XL,XU,D0,D1 , Cl ,DU,DL 

DX=(PMAX-PHIN)/10.0 
XL=PMIN-D X 
XU=PMAX+DX 
DO=ABS(PMIN) 

D1=ABS(PMAX) 

IF(D1  .LT.  DO)  D1=D0 
NN=21 

DO  100  JJ=1 ,40 
NNsNN-1 
C1=10.0»»NN 

IF(D1  .GE.  Cl)  GOTO  200 
100  CONTINUE 
200  CONTINUE 
D1=PMAX/C1 
II=IFIX(D1) 

DU=FLOAT(II)*C1 
IF (DU  .LT.  PMAX )  11=11+1 

DU=FLOAT(II)*C1 
IF (DU  .GT.  XU)  THEN 
C1=C1/10.0 
GOTO  200 
END  IF 

D1=PMIN/C1 
II=IFIX(D1 ) 

DL=FLOAT(II)»C1 
IF(DL  .GT.  PMIN)  11=11-1 
DL=FL0AT(II)«C1 
IF(DL  .LT.  XL)  THEN 
C1=C1/10.0 
GOTO  200 
END  IF 
PMINrDL 
PMAX=DU 

RETURN 

END 
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SUBROUTINE  BORDER  ( JK , IK , XLAB , TL AB , XT , YT , BUP ) 

Subroutine  to  provide  borders  for  printer  plotting 

CHARACTER* 1  XT( 8 ) , IT ( 8 ) , BUP ( 1 28 ) , TEMP (10) 

INTEGER  I,J,JK,IK,ITEMP,NUMLAB 
REAL  XLAB(1 1 ),TLAB(1 1 ) 

PARAMETER  (no-0) 

PARAMETER  (yes=1) 

NUMLABsno 

IF  (IK  .BQ.  0  .OR.  MODUK.5)  .EQ.  0)  THEN 

NUMLAB=yes  I  Numeric  label  required 

JK=JK  1 
END  IF 

IF (IK  .EQ.  0  .OR.  IK  .EQ.  50)  THEN  I  Top  or  bottom  border 

ENCODE  (10, 900, TEMP)  YLAB(12-JK) 

DO  100  1=1,10 

BUF(I+12)=TEMP(I) 

100  CONTINUE 

DO  300  1=2,11 

DO  200  J=6, 14 

ITEMP=I*10  ♦  J 

IF(BUF(ITEMP)  .EQ.  '  ’)  BUF(ITEMP)= 

200  CONTINUE 

ITEMP=I*10  +  5 

IF(BUF(ITEMP)  .BQ.  •  ')  BUFdTEMPjs'I' 

300  CONTINUE 

NUMLAB=no 

ELSE 

IF (IK  .BQ.  25)  THEN  I  Title  for  vertical  axis 

DO  400  1=1,8 

BUF(I+2)=YT(I) 

400  CONTINUE 

END  IF 
END  IF 

IF(NUMLAB  .EQ.  yes)  THEN  1  Numeric  labels  and 

ENCODE  (10, 900, TEMP)  YLAB(12-JK)  1  associated  symbols,  but 
DO  500  1=1,10  I  not  for  top  &  bot.  borders 

BUP(I+12)=TEMP(I) 

500  CONTINUE 

IF(BUF(25)  .EQ.  »  ')  BUF(25)  ='■*•’ 

IF(BUF(26)  .EQ.  '  »)  BUF(26)  =  •-• 

IF(BUF( 124 )  .EQ.  '  ')  BUF(124)=*-' 

IF(BUF(125)  .EQ.  1  »)  BUF( 125)= 

ELSE 

IF(BUF(25)  .BQ.  '  •)  BUF(25)  ='I' 

IF(BUF( 125)  .BQ.  •  ')  BUF(125)='I' 

END  IF 

WRITE (6, 904)  (BUF(I) ,1=1 ,125)  I  Print  contents  of  buffer 

IF (IK  .EQ.  50)  THEN 

CALL  BZER0  (BUF.128)  I  Get  labels  for  bottom  border 

DO  600  1=1,11 

ENCODE  (9, 902, TEMP)  XLAB(I) 

DO  600  J=1 ,9 

ITBMP=I*10  ♦  9  ♦  J 
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C 

900  F0RMAT(F10.4 ) 
902  FORMAT  (F9.4) 
904  FORMAT ( 128A1 ) 
906  FORMAT(/) 
RETURN 
END 
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SUBROUTINE  SETCRD  (IK, TEMPI ,ICXY) 

Subroutine  to  compute  rounded  (Integer)  data  point  coordinates 

INTEGER  IK , ICZY ( 1 ) 

REAL  TEMPI , TBMP2 
C 

IF (TEMPI  .LT.  0.0)  TEMPIsO.O 
TEMP2=AINT( TEMPI ) 

IF ( TEMPI  -  TBMP2  .GT.  0.5)  TEMP2=TEHP2  +  1.0 
ICZT ( IK )=IFIX (TEMP2 ) 

C 

RETURN 

END 
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SUBROUTINE  BZERO  (B,N) 


Subroutine  to  initialize  a  character  array 

CHARACTER* 1  B(N) 

INTEGER  I , N 
C 

DO  100  1=1, N 
B(I )= '  ' 

100  CONTINUE 
RETURN 
END 
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EVAL 


C 
C 

C  Program  to  predict  homogeneous  test  results  for  material  models. 

C  Prepared  by  L..R.  Herrmann  at  the  University  of  California, 

C  Davis  campus.  Fortran  77  version,  revised  1983. 

C 

INTEGER  I , J , K , NJ , K7 , ITNO , NUM , INC , NINC ,  IPRNT1  ,  IPRNT2 ,  UMAX  ,  KIND, 

*  LARGE, LOCI T, JUNK, If LAG, I PR I NT , I TITLE (40 ) ,IC0D(7 ) ,IPL0T( 10) 
REAL  PR0P( 19), STOR ( 7 ) , SIGB( 6 ) , DSIG(6 ) , EPB(6 ) , DEP(6 ) , V (6 ) , DV ( 6 ) , 

*  RP (6 ) , R (7 ) » DB(6 , 6 ) , DE(6 , 6 ) , S(7 , 7  ) , PINC ( 1 50 , 8 } , CONFIN , 

*  ERMAX  f  CONFL , SNORM 1 f SN0RM2 , ENORM 1 , EU0HM2 , DLTA U , U , D , GAM , 

»  TEMP, TEMPI , TEMPI T 

PARAMETER  (yes=1) 

PARAMETER  (no=0) 

C 

1  CONTINUE 

READ (5,800, END=999 )  (ITITLE(I ) ,1=1 ,40) 

WRITE ( 6 , 900 )  (ITITLE(I ) ,1=1,40) 

READ( 5 , 802 )  ST0P{ 1 ) , STOR( 7 ) , CONFIN 
C 

C  Read  material  properties 

C 

CALL  RPROP  (PROP) 


C 

C 


c 


c 

c 

c 


Read  iteratcn  and  analysis  information 

READ (5,801)  ITMAX , ERMAX , CONFL , KIND , LARGE , LOCIT 
IF(ITMAX  .GT.  20)  ITMAX=20 
IF (CONFL  .LE.  0.0)  CONFL=0.3 
IF ( LOCIT  .EQ.  0)  LOCIT=1 


R£AD(5, 804)  IPRNT1 , IPRNT2  l  Read  flags  for  output  control 

READ (5, 804)  (JPLOTU ) , 1=1 , 10)  !  Read  plotting  instructions 


WRITE (6 ♦ 902 )  STOR (7 ) . CONFIN , STOR ( 1 ) 

STOR(1 )=ST0R(1 )f3«0 

IF (KIND  .EQ.  0)  WRITE (6 , 906 ) 

IF (KIND  .EQ.  1)  WRITE(6,907) 

IF(LARGE  .EQ.  0)  WRITE(6,908) 

IF(LARGE  .EQ.  1)  WRITE(6,909) 

WRITE(6,90 ! )  ITMAX, LOCIT, ERMAX, CONFL 


Initia  1  izat  ion 


NUM=0 

IFLAG=no 

IPRINT=no 

ICOD(7)=0  !  Fictitious  increment  type  specification 

TEMPI T=  0 . 0 

SN0RM1 =0, 0 

ENORM1 =0.0 

DLTAU  =  0 . 0 

U-  0 . 0 

DO  100  Isl f 6 

SIGB(I )  =  0 . 0 
DSIG( I )  =  0 . 0 
EPB( I )=  0.0 
DEP(I )=  0.0 

100  CONTINUE  59 
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DO  150  1=1,3 

SIGB(I)=C0NFIN 
150  CONTINUB 

Loading  segment  loop 

IF(IPRNT1  .BQ.  0)  WRITE (6, 903) 

200  CONTINUE 

RE AD (5, 803)  (ICOD(I) ,T(1) ,1=1 ,6) ,NINC,D 
IF(IC0D(1)  .GT.  1)  GOTO  750 
IF(D  .BQ.  0.0)  D=1.0 

Determine  first  Increments 

TEMPI = 1 . 0/ ( FLOAT ( NINC ) ) 

IF(D  .NB.  1.0)  TO!P1  =  ( 1.0  -  D)/(1.0  -  D»»NINC) 
DO  250  1=1,6 

TEMP=V ( I )  -  SIGB(I) 

IF(ICOD(I)  .BQ.  1)  TEMP=V(I)  -  EPB(I) 

DV ( I ) =TEMP*TEMP 1 
250  CONTINUE 


Change  sign  of  strain  estimate  at  beginning  of  new  segment 
in  case  of  unstable  behavior  at  the  end  of  the  previous  one. 

DLTAU=DLTAU#0. 01 
DO  300  1=1,6 

DSIG(I)=DSIG(I)».01 
DEPd)=DEPd)*-.01 
300  CONTINUB 

Increment  loop 

DO  700  INC=1 , NINC 


Iteration  loop  (successive  approximation) 

DO  600  ITN0= 1 , 1TMAX 

Get  Incremental  properties 

NJ=NUM  +  1 
K7=ITN0 

CALL  CLAY  (3,NJ,E7,BRMAX,PROP,STOR,SIGB,EPB,DSIG,DEP, 
•  DB, DE,U,DLTAU, GAM, KIND, LARGE, L0CIT, 0.5) 

Form  and  modify  stiffness 

DO  350  1=1,3 
S(7,I)=GAM 
S(7,I+3)=0.0 
S(I,7)=1.0 
S(I+3,7)*0.0 
350  CONTINUB 

S(7,7)=-1.0 
R ( 7 )*0.0 
DO  450  1=1,6 

DO  400  J=1 ,6 

S(I,J)=0.5»(DB(I,J)  ♦  DE(I, J) ) 

400  CONTINUB 
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R(I)=DV(I) 

IF(ICOD(I)  .EQ.  1)  THEN 
DO  420  K=1 ,7 

IF(ICOD(K)  .EQ.  0) 

•  R(K)=R(K)  -  S(K, I)»DV(I) 

S(I,K)=0.0 

S(K,I)=0.0 

420  CONTINUE 

S(I,I)=1.0 
END  IF 

450  CONTINUE 


Solve  for  strain  increment 
CALL  SOLVE  (KIND,IC0D,S,R) 


Calculate  total  stress  increments 

DO  550  1=1,6 
TEMP=0. 0 

IF(I  .LT.  4)  TEMP=R(7 ) 

DO  500  J= 1 , 6 

TEMP=0. 5*(DB(I , J )  ♦  DE(I,J))*R(J)  +  TEMP 
500  CONTINUE 

RP(I)=TEMP 

550  CONTINUE 

Compute  error  norms 

CALL  ERNORM  (IFLAG.IPRINT, ITN0,NUM,KIND,ERMAX,C0NFL, 

*  SN0RM 1 , SNORM2 , EN0RM1 , EN0RM2 , TEMP 1 T , RP , 

*  DSIG , R , DEP , DLTAU ) 

IF(IFLAG  .EQ.  yes)  GOTO  640 

600  CONTINUE 

IPRINT=yes 

CALL  ERNORM  (IFLAGdPRINT, ITNO.NUM, KIND, ERMAX, CONFL, SNORM1 , 

*  SN0RM2 , ENORM 1 , EN0RM2 , TEMP 1 T , RP , DSIG , R , DEP , DLTAU ) 

Upon  failure  to  converge,  eat  up  remaining  input  data 

620  READ (5, 803)  JUNK 

IF( JUNK  .GT.  1)  GOTO  750 
GOTO  620 

Update  total  values 

640  NUM=NUM  +  1 

TEMP1T=0.0 
DO  660  1=1,6 

DV(I)=DV(I)»D 
DSIG(I )=RP(I ) 

DEP(I)  =R(I) 

SIGB(I )=SIGB( I )  +  DSIG(I ) 

EPB(I)  =EPB( I )  +  DEP(I) 

TEMPI T  =TEMP IT  ♦  ABS(SIGB(I ) )»0. 1 
660  CONTINUE 

IF(KIND  .EQ.  0)  DLTAU=R(7 ) 

U=U  +  DLTAU 

Store  incremental  values  for  future  plotting 


PINC(NUM,  1  )  =  (SIGB(3)  -  SIGB(1)) 

PINC(NUM,  2)=(SIGB(1 )  ♦  SIGB(2)  ♦  SIGB(3))/j-0  -  U 
PXNC(NUM  3 ) r  U 

PINC(NUm’,  4)=(EPB( 1 )  +  EPB(2)  +  EPB(3»§100,0 
PINC(NUM,  5)=(EPB(3)  -  EPB(1))/(1.50)«100.0 
PINC (NUM,  6)  =  BPB(3)§100.0 

PINC(NUM,  7)=(SIGB(3)  -  SIGB(1 ) )/(STOR( 1 )/3.0) 

PINC(NDH,  8)  =  (SIGB(-1 )  ♦  SIGB(2)  ♦  SIGB(3)  -  3.0»U)/STOR(1 ) 


Print  incremental  values  of  stresses  and  strains  if  desired 
IF(IPRNT1  .EQ.  0) 

•  WRITE (6 ,904)  NUM, (EPB(I) ,1=1 ,6) , (SIGB(I) ,1=1 ,6 ) .U.ITNO 
700  CONTI NOE 

WRITE (6, 905) 

GOTO  200 
750  CONTINUE 

Print  computed  incremental  parameters  if  desired 

IF(IPRNT2  .EQ.  0)  THEN 
WRITE(6,910) 

DO  780  1=1, NUM 

WRITE(6 ,91 1 )  I,PINC(I ,6 ) , (PINC(I, J) , J=1 ,5) 

780  CONTINUE 
END  IF 

IF(IPLOT(  1)  .EQ.  yes) 

•  CALL  JPL0T( NUM, PINC, 2,1, *  P  Q  ') 

IF(IPLOT(  2)  .EQ.  yes) 

•  CALL  JPLOT(NUM,PINC,6, 1 , '  EPS-1  Q  ') 

IF(IPL0T(  3)  -EQ.  yes) 

•  CALL  JPLOT(NUM,PINC,6,3, *  EPS-1  V  0  ') 

IF (I PLOT (  4)  .EQ.  yes) 

•  CALL  JPLOT(NUM,PINC,6,4, '  EPS-1  E-VOL  ’) 

IF ( I PLOT (  5)  .EQ.  yes) 

•  CALL  JPL0T( NUM, PINC, 8,7, '  P/PO  Q/PO  ') 

IF( IPL0T(  6)  .EQ.  yes) 

•  CALL  JPLOT(NUM,PINC,6,7, '  EPS-1  Q/PO  •) 

GOTO  1 

999  CALL  EXIT 

Format  statements 

800  FORMAT ( 4 0A2) 

801  F0RMAT(I5,2E10.3,3I5) 

802  F0RMAT(3E10.3) 

803  FORMAT(6(I1,B9.1),I5,E10.3) 

804  FORMAT( 1015) 

900  FORMAT ( 1H1 ,/,4X,40A2,4(/) ) 

901  FORMAT (/,1 OX, 'ITERATION  AND  CONVERGENCE  PARAMETERS:',/, 

1  24X , ’ ITMAI= ' , 13 , / • 24X , *LOCIT= • ,13,/, 

2  24X, 'BRMAXs* ,F6.3»/ »24X, 'C0NFL=' ,F6.3,/) 

902  FORMAT (/,29X,* INITIAL  VOID  RATIO  =\F7.4,/, 

1  21X,' INITIAL  CONFINING  PRESSURE  =' , 1PE12.3,/, 

2  10X, 'INITIAL  PRECONSOLIDATION  PRESSURE,  PO  =',1PE12.3,//) 

903  FORMAT ( 1 H 1 , 3X , ' N ' , 4X , ' BPS-X ' , 4X , ' BPS-Y • , 4X , ' BPS-Z • , 3X , ' GAM-XY • , 

1  3X, 'GAM-XZ' ,3X, *GAM-YZ*,5X, 'SIG-X'.SX, ’SIG-Y',5X, 'SIG-Z’, 

2  4X , * TAU-XY ' , 4X , ' TAU-XZ ' , 4X , ' TAU- YZ ' , 6X , ’ U ' , 3X , ' IT  #',/) 


904  FORMAT(1X,I3,1P6E9. 1.7E10.2, IX, 13) 

905  FORMAT(/) 

906  FORMAT (5X, ’•••**  REFORMULATED  NEARLY- INCOMPRESSIBLE 

1  'ANALYSIS  ••*•*'/) 

907  FORMAT (5X, '•*•••  NON-REFORMULATED  NEARLY-INCOMPRESSIBLE 

1  'ANALYSIS  *••••'/) 

908  FORMAT ( 5X , ' •**••  ENGINEERING  STRESSES  AND  STRAINS 

1  'ASSUMED  •••••'/) 

909  F0RMAT(5X, '•••••  TRUE  STRESSES  AND  NATURAL  STRAINS 

1  ’ASSUMED  **•••'/) 

910  FORMAT ( 1H1 , //, 9X, ' N’ ,5X, 'EPS-1 ',9X, 'Q' ,9X, 'P' ,9X, ’ U ' , 

1  5X,  'E-VOL'  ,5X,  'E-DEV' /9X, '-'  ,5X, ' - '  ,5X, ' - ' , 

2  5X, ' - '  ,5X,  ' - \5X, ' - ’  ,5X, ' - '/) 

911  FORMAT (5X,I5,6F10.2) 

END 


63 


ooo  oooooo 


SUBROUTINE  SOLVE  (EIND,ICOD,S,R) 

Subroutine  to  solve  tbe  7x7  stiffness  matrix  for  EVAL 

INTEGER  I,II,IL,J,E,NfKIND,IC0D(7) 

REAL  TEM?,R(7),S(7,7) 

Reduction  of  stiffness  matrix  and  r.h.s.  vector 

N=7-KIND 
DO  400  1*1, N 

IF(IC0D(I)  .EQ.  0)  THEN  !  Avoid  work  for  trivial  rows 

TEMP* 1 .0/S(I,I) 

R(I)sR(I)*TEHP 
DO  1 00  J=I,N 

S(I,J)=S(I,J)«TEMP 
100  CONTINUE 

IF (I  .NE.  N)  THEN 
II=It-1 

DO  300  J=II,N 

IF(IC0D(I)  .EQ.  0)  THEN 
TEMP=-S(J,I) 

DO  200  K=I,N 

S(J,K)=S(J,K)  +  TEHP*S (I , K ) 

200  CONTINUE 

R(J)=R(J)  +  TEMPER (I ) 

END  IF 

300  CONTINUE 

END  IP 
END  IF 

400  CONTINUE 

Back  substitution 

I=N 

DO  500  11=2, N 
1=1-1 
IL=I+1 

DO  500  J=IL,N 

R(I)=R(I)  -  S(I,J)»R(J) 

500  CONTINUE 
RETURN 
END 
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SUBROUTINE  ERNORM  ( IFLAG , IPRINT, ITNO , NUM, KIND , ERMAX , CONFL , SN0RM1 , 

•  SN0RM2 , ENORM 1 , EN0RM2 , TEMP 1 T , RP , DSI G , B , DEP , DLTA  U ) 

Subroutine  to  check  for  convergence  of  the  incremental  solution 
I NTEGER  I , IFL AG , I PRI NT , I TNO , NUM , K I ND 

REAL  ERSIG , EREPS , ERMAX , RP ( 6 ) , DSIG ( 6 ) , R ( 7 ) , DEP ( 6 ) , DLTAU , TEMP 1 , 

*  TEMP2 , TEMP 1 T , CONFS , CONFE , CONFL , SN0RM1 , SNORM2 , ENORM 1 , EN0RM2 
PARAMETER  (yes=1) 

PARAMETER  (no=0) 

ERSIG=0.0  !  Compute  error  norms 

EREPS=0. 0 
TEMPUO.O 
TEMP2=0. 0 
DO  100  1=1,6 

ERSIG=ERSIG  +  ABS(RP(I)-DSIG(I) ) 

EREPSrEREPS  +  ABS(R(I )-DEP(I) ) 

TEMPI =TEMP1  +  ABS(RP(I)) 

TEMP2=TEMP2  +  ABS(R(I)) 

100  CONTINUE 

IF (TEMPI  .LT.  TEMP1T)  THEN 
ERSIG=ERSIG/TEMP1T 

ELSE 

ERSIG=ERSIG/TEMP1 
END  IF 

EREPS=EREPS/TEMP2 

Upon  failure  to  converge,  print  pertinent  information  about  errors 

IF(IPRINT  .EQ.  yes)  THEN 

WRITE (6,900)  NUM+ 1 , ERSIG , EREPS , ERMAX 
RETURN 
END  IF 
C 

C  Check  norms  against  error  tolerance  to  determine  convergence 

C 

IF (ERSIG  .GE.  ERMAX  .OR.  EREPS  .GE.  ERMAX)  THEN 
C 

C0NFS=1.0  !  Apply  Aitken’s  convergence  acceleration 

CONFE= 1 . 0 

IF (ITNO  .NE.  1  .AND.  (-1)»»ITNO  .LT.  0)  THEN 
CALL  ACCEL  (SN0RM2.SN0RM1 , TEMPI , CONFS, CONFL) 

CALL  ACCEL  ( EN0RM2, ENORM 1 ,TEMP2, CONFE, CONFL) 

END  IF 
TEMP1=0.0 
TEMP2=0. 0 
DO  200  1=1,6 

DSIG(I )=RP(I )*C0NFS  +  ( 1 .0-CONFS)«DSIG(I ) 

DEP(I)  =R(I )*C0NFE  +  ( 1 . 0-CONFE)*DEP(I ) 

TEMP1=TEMP1  ABS(DSIG(I)) 

TEMP2=TEMP2  +  ABS(DEP(I)) 

200  CONTINUE 

IFOCIND  .EQ.  0)  DLTAU=R(7)*CONFS  +  ( 1 .0-C0NFS)«DLTAU 
SN0RM2=SN0RM1 
EN0RM2=EN0RM1 
SN0RM1=TEMP1 
EN0RMUTEMP2 
IFLAG=no 

ELSE 
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IPLAQsyes 
END  IF 
C 

900  FORMAT (///,3X, ’•••*•  CONVERGENCE  DID  NOT  OCCUR  FOR  INCREMENT  ',13, 

1  t  ••••«.  f/j0x,  'ERSIGs’.IPEIO.S.SXf'EREPSs'.IPEIO.S, 

2  5X, 'ERMAXs* , 1  PEI 0.3*///) 

RETURN 

END 
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SUBROUTINE  ACCEL  (X2,X1 ,Xf C,XL) 

C 

C  This  subroutine  calculates  the  Aitken's  convergence  factor 
C 

REAL  X,X1 ,X2,XL,C,TEMP 
C=1.0 

TEMP=-X2  -*•  2.0*X1  -  X 
IF (TEMP  .NE.  0.0)  THEN 
C=(X1  -  X2)/TEMP 
IF(C  .LT.  XL)  C=XL 

IF(C  .GT.  1.0/XL)  C=1 . 0/XL 

END  IF 
RETURN 
END 
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SUBROUTINE  RPROP  (PROP) 

This  subroutine  reads  in  and  modifies  the  parameters  required 
by  the  bounding  surface  plasticity  model  for  cohesive  soils. 

REAL  PR0P(19) 

READ ( 5 , 800 )  (PROP(I),ISl,3)»(PROP(I),l=0,1l),PROP(7),PROP(5), 

1  PROP(8),PROP(6),  PROP (17),  PR0P(16),  PR0P(19), 

2  PROP (4 ) , PR0P( 18) , (PROP(I) , 1=12,15) 

WRITE(6,900)  PR0P(1),  PR0P(3),  PR0P(2),  PR0P(4) 

IF(PR0P(5)  .LT.  0.5)  WRITE(6,901)  PR0P(5) ,PR0P(7) 

IF(PR0P(5)  .GE.  0.5)  WRITE(6,902)  PR0P(5),PR0P(7) 

WRITE(6,903)  PR0P( 17) ,PR0P( 19) ,PR0P( 16) ,PR0P( 18) 

WRITE (6, 904 )  PROP(9),  PR0P( 12) ,PR0P( 10) ,PR0P( 13),PR0P( 11 ) 

Check  the  magnitude  of  the  shape  parameter  * T " 

CALL  TCHECK  (PROP) 

WRITE ( 6 , 905 )  PROP (15), PROP ( 1 4 ) , PROP ( 8 ) 

IF(PR0P(6)  .EQ.  0.0)  WRITE(6,906) 

IF(PR0P(6)  .NE.  0.0)  WRITE (6, 907)  PR0P(6) 


Convert  parameters  from  triaxial  to  invariant  stress  space 

PR0P(3)  =  PROP ( 3 ) /SQRT (27.0) 

PR0P(7)  =PROP(7)*3.0 
RETURN 


Format  statements  for  RPROP 
800  FORMAT(8E10.3) 

900  FORMAT (1  OX, 'CLASSICAL  CLAY  MATERIAL  CONSTANTS:',/, 

1  15X, 'LAMBDA  = • ,F7.4 , 17X, 'MC  =',F7.4,/, 

2  15X, 'KAPPA  =' ,F7.4, 14X, 'ME/MC  =',F7.4,/) 

901  F0RMAT(15X, ’POISSON»S  RATIO  = ' ,F7.4, 8X, 'PL  = ' , 1PE10.3,//) 

902  F0RMAT(15X, 'SHEAR  MODULUS  = ' , 1PB10. 3.7X, 'PL  = • , 1PE10.3,//) 

903  FORMAT (1 OX, ’HARDENING  PARAMETERS:',/, 

1  1 5X , ’ SM  =',F7.4,21X,’H2  =’,F7.4,/, 

2  1 5X , ' HC  = ' , F7 . 4 , 1 8l , ' HE/HC  =',F7.4,//) 

904  FORMAT(1 OX, 'PARAMETERS  DESCRIBING  SHAPE  OP  BOUNDING  SURFACE:',/, 

1  151, ’RC  =’ ,F7.4, 18X, 'RE/RC  =',F7.4,/, 

2  1 51 , ' AC  = ' , F7 . 4 , 1 81 , ' AE/AC  =’,F7.4,/, 

3  16X,'T  =’,F7.4,//) 

905  FORMAT (17X, 'ELASTIC  NUCLEUS  PARAMETER,  S  =',F7.4,/, 

1  15X, 'PROJECTION  CENTER  PARAMETER,  C  =’,F7.4,//, 

2  251, 'ATMOSPHERIC  PRESSURE  = * , 1  PE 10. 3,/) 

906  FORMAT (//,25X, '•••»•  DRAINED  CONDITIONS  •••••',/) 

907  F0RMAT(//,5X, *•••••  UNDRAINED  CONDITIONS  -  THE  COMBINED  ', 

1  'SKELETON  AND  WATER  BULK  MODULUS  =\1PE10.3,'  •••••',/) 

END 
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SUBROUTINE  TCHECK  (PROP) 

This  subroutine  checks  the  value  of  the  bounding  surface  shape 
parameter  "T"  and  adjusts  this  value  if  it  exceeds  the  theoretical  max. 
Original  version  written  by  J.S.  De  Natale. 

REAL  TEMP 1 , TEMP2 , TEMPR , PROP (19) 

YFUN(TT)=(1 .0  ♦  TT)*SQRT(1 .0  ♦  TT«TT)  -  (1.0  +  TT«TT) 

Check  against  theoretical  limit  in  compression 
TEMPI =PR0P(11) 

TEMP2= PROP ( 9 ) •PROP (10) »SQRT (27.0)/ PROP ( 3 ) 

TEMP2= YFUN ( TEMP2 ) 

TEMP2=TEMP2/2 . 0/PR0P (9 ) 

IF ( PROP (1 1 )  .GT.  TEMP2)  PR0P( 1 1 )=TEMP2 


Check  against  theoretical  limit  in  extension 
TEMPR=PR0P ( 9)*PR0P (12) 

TEMP2=TEMPR»PR0P ( 1 0 ) *PR0P (13) *SQRT (27.0 ) /PROP( 3 )/PR0P ( 4 ) 

TEMP2=  YFUN ( TEMP2 ) 

TEMP2=TEMP2/2 . 0/TEMPR 

IF ( PROP ( 11)  .GT.  TEMP2 )  PROP( 1 1 )=TEMP2 
IF(PR0P( 11)  .NE.  TEMPI)  WRITE(6,100)  PR0P(11) 

100  F0RMAT(5X, *•••*•  THE  USER-SPECIFIED  VALUE  OF  T  EXCEEDS  THE  MAX', 

1  '  PERMISSIBLE  VALUE  ,/,5X, '•••••  T  HAS  THUS', 

2  '  BEEN  AUTOMATICALLY  RESET  TO  T  =',F7.4,8X,'  •••••',//) 

RETURN 

END 
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